Moving NRQCD for heavy-to-light form factors on the lattice 
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We formulate Non-Relativistic Quantum Chromodynamics (NRQCD) on a lattice which is boosted 
relative to the usual discretization frame. Moving NRQCD (mNRQCD) allows us to treat the 
momentum for the heavy quark arising from the frame choice exactly. We derive mNRQCD through 
0{l/m^ ,vfg,i), as accurate as the NRQCD action in present use, both in the continuum and on the 
lattice with 0{a*) improvements. We have carried out extensive tests of the formalism through 
calculations of two-point correlators for both heavy- heavy (bottomonium) and heavy-light (Bs) 
mesons in 2+1 flavor lattice QCD and obtained nonperturbative determinations of energy shift 
and external momentum renormalization. Comparison to perturbation theory at 0{as) is also 
made. The results demonstrate the effectiveness of mNRQCD. In particular we show that the decay 
constants of heavy-light and heavy-heavy mesons can be calculated with small systematic errors up 
to much larger momenta than with standard NRQCD. 
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I. INTRODUCTION 

The Cabibbo-Kobayashi-Maskawa (GKM) matrix is 
the focus of intense study; an inconsistency between in- 
dependent determinations of GKM matrix elements from 
different physical processes would be evidence for new 
physics beyond the Standard Model. While experimen- 
tal measurements of exclusive semileptonic decays have 
reached good precision and will be improved further by 
LHGb, determinations of CKM matrix elements from the 
decay rates are complicated by the need for precise theo- 
retical calculations in nonperturbative quantum chromo- 
dynamics (QCD). Lattice QCD provides a first-principles 
approach to these calculations and it is important to re- 
duce systematic and statistical errors as far as possible. 

For example, the decay B iriv [l|, d, Q can be 
used to determine the CKM matrix element Vub while 
the rare decays B K*-f, K(*H+1- 0, d, i, 0, pro- 
vide excellent opportunities to study contributions from 
new physics, as the flavor-changing neutral current 6 — > s 
is loop-suppressed in the Standard Model. In both cases, 
a nonperturbative calculation of the hadronic form fac- 
tors is required. 

These form factors are a function of the momentum 
transfer squared, g^, where q = pb ~ Pf is the difference 
between the four-momenta of the B meson and the meson 



in the final state. If this meson is light compared to 
the B meson, the recoil momentum at small values of 
can be very large. Unfortunately, current lattice QCD 
calculations of these form factors work well only for low 
recoil momenta, i.e. large g2 d, [H [H [il], while for 
B K*j one has q^ = and experimental data for 
B 7r£i/ covers the full q-^ range 

By computing at just one or a few points with large 
g^, one might be able to reduce the error on \Vub\ from 
B -kIv, where the shape of the form factor is now be- 
ing measured precisely by experiment [J, [3, 0]- However, 
the form factors governing the rare 6 — *■ s decays are 
not well-determined and must be computed using lattice 
QCD. Given the propensity for models of new physics 
to introduce new sources of flavor-changing neutral cur- 
rents, it is desirable to have new tools to reduce the errors 
on the Standard Model calculations of differential cross 
sections for rare decays. 

In this paper we present a technique for extending lat- 
tice QCD calculations of the decays of mesons contain- 
ing one heavy quark to lower q^ values than has hitherto 
been possible by reducing the discretization errors owing 
to the large recoil of the final state meson. 

The formalism that we describe and put to the 
test in subsequent sections is a generalization of Non- 
Relativistic QCD (NRQCD) The NRQCD for- 
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malism, which has already had considerable success in 
the study of heavy-quark systems, relies on the fact 
that fluctuations in the heavy quark momentum within 
a heavy meson are small compared with the mass of the 
meson itself. The Lagrangian of NRQCD is expressed 
as a sum over operators whose importance is governed 
by power-counting rules; in dimensionless units the op- 
erators are ordered in powers of g, Wroi for heavy-heavy 
mesons and in powers of asjAgcD/'Ti for heavy-light 
mesons, where = j (47r) is the strong couphng con- 
stant, Wrci ~ IpI/'ti is the relative internal velocity of 
the heavy quarks and to is the heavy quark mass.^ For 
NRQCD, the heavy meson is usually taken to be at rest 
in the lattice frame. This is appropriate for calculations 
of the mass spectrum of heavy-light and heavy-heavy 
mesons and for zero-recoil or low-recoil decays. How- 
ever, for the heavy-to-light decays of the _B-meson cited 
above, outside the low recoil region the momentum of 
the light meson in the final state becomes comparable 
to the inverse lattice spacing. Consequently the calcu- 
lation is sensitive to lattice artifacts which lead to large 
systematic errors. 

It is therefore better to give the B meson a non-zero 
momentum in the opposite direction, thereby reducing 
the final meson's momentum at a given . To substan- 
tially reduce the momentum of the final meson, the mo- 
mentum of the B meson has to be very large, so that 
NRQCD would no longer be able to describe the h quark 
inside it due to relativistic and lattice spacing errors. 
However, we note that fluctuations of momentum of the 
heavy quark inside the B meson are much smaller than 
the momentum of the meson itself. Therefore, to reduce 
errors, instead of discretizing the momentum of the h 
quark itself, we choose to discretize its fluctuations in- 
side the moving B meson. The formalism which achieves 
this goes by the name of moving NRQCD (mNRQCD) in 
which the expansion is about the state where the heavy 
quark is moving with a velocity w, the frame velocity; this 
formalism was introduced briefly in [Tsj : Earlier, related 
approaches were proposed in [3, . 

The remainder of the paper is structured as follows. In 
Section |TT] we discuss the choice of the optimal reference 
frame for the lattice calculations. We give an explicit 
derivation of the continuum mNRQCD action in Sec- 
tion [Hll We explain how the theory is discretized in Sec- 
tion [IVl In Section |V] we develop perturbative methods 
for mNRQCD and explain how to derive the renormaliza- 
tion of parameters due to radiative corrections. We give 
1-loop results for the heavy quark renormalization con- 
stants. The construction of decay currents is discussed 
in section IV Dl Then, in section IVII we present the re- 
sults of nonperturbative calculations based on two-point 



^ These rules are frequently referred to as NRQCD and HQET 
power counting schemes, respectively. Note that the choice of 
NRQCD as a lattice action is compatible with both schemes. 
See Sec. HITB] below. 



correlators for heavy-heavy and heavy-light mesons in 
mNRQCD. These include the spectrum, renormalization 
constants and decay constants for various values of the 
frame velocity v. 

The perturbative and nonperturbative renormalization 
constants are compared in Section IVIII Wc summarize 
and discuss our results in Section IVlIII 

In the Appendices we specify some notation (Ap- 
pendix 1^ , describe the removal of time derivatives in 
the ©(A^cd/to^) mNRQCD Hamihonian ( Appendix E]), 
give explicit expressions for the lattice derivative opera- 
tors (Appendix [C]) and tadpole improvement corrections 
(Appendix [D| and present further perturbative results 
for a set of simpler actions (Appendix [E|) . We comment 
on the poles of the Symanzik-improved gluon action in 
Appendix iFl 

Preliminary versions of this work have been presented 
in Refs. [H, [3, S [HI, S HI • 

II. MINIMIZING ERRORS 

We start by parametrizing the 4-momcntum of the h 
quark as 

p = TO u -f A: 

where to is the mass of the 6-quark, and m a 4- 
velocity. In traditional (non-moving) NRQCD one has 
u — (1,0,0,0), and a non-relativistic expansion in the 
residual 3-momcntum fc is performed. In other words, the 
heavy-quark mass term is removed from the Lagrangian. 
Thus, the 3-momentum p, which is equal to fc in this 
case, has to be small to prevent large relativistic errors 
as well as discretization errors on the lattice. 

In moving NRQCD, we generalize this to other frames 
of reference, removing the momentum to u with an arbi- 
trary 4-velocity u from the Lagrangian, and again per- 
forming a non-relativistic expansion in the residual 3- 
momentum fc. The relativistic energy of the heavy quark 
is -E = y^j^ + TO? = ^ [mu -f fe)2 -|- w? . Taylor expan- 
sion for small |fc| gives 

. -(v fc)2 
= 7TO + t) • fe H — H 

27TO 

where we write u = (uq, w) — (7, ^v) with the 3- velocity 
V and 7 = (1 — v^)^^!"^ . Discarding the constant term 
7TO, we expect that the 0(1/to) "kinetic" part of the 
continuum mNRQCD Hamiltonian in momentum space 
will be given by 

Ha=v-k + ^ '-. 

27TO 

Of course, the size of k and the associated relativistic 
and discretization errors depend on the choice of u. The 
standard choice is m = ps/Mg, the 4-velocity of the B 
meson. Then, the residual momentum k is small com- 
pared to pg and the non-relativistic expansion in k is 
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a good approximation even for B mesons at moderately 
high velocities. 

1. Discretization errors 

One of the main applications of the niNRQCD ap- 
proach is to the heavy-to-light weak decay of a i?-meson 
to a final state including a light meson. As discussed in 
the introduction the size of the discretization errors in a 
lattice calculation depends on the momentum of the fi- 
nal state meson; states with spatial momenta comparable 
to the inverse lattice spacing can be affected by lattice 
artifacts. Nevertheless, one wishes to compute matrix el- 
ements over the whole physical kinematic range, includ- 
ing the large recoil regime where the final state has large 
momentum relative to the B meson. With mNRQCD 
we attempt to reduce discretization errors by choosing 
a non-zero frame velocity w, so that the final state me- 
son can have moderate spatial momentum in the lattice 
frame, even as we explore large recoil kinematics. 

If the B meson is at rest, the residual momentum k 
has a distribution with width of the order Aqcd and the 
residual energy has a distribution with width of the order 
^QCd/(2™) Aqcd- Note that the momentum p^^^^ of 
the light quark in the B meson (the "spectator quark") 
is of the same order by momentum conservation. 

For a B meson moving with velocity v, the momentum 
distribution is boosted to approximately 7 Aqcd- Let us 
now consider a decay B ^ F where F denotes the light 
meson in the final state and the 4-momenta are 

PB = {-fMB, iMbv), 

PF = {'jMj, + \pp\^, pp) 

where v is antiparallel to Pp. For a given value of 
9^ = [pb — PfY , we shall determine the optimal veloc- 
ity of the B meson which minimizes discretization errors. 
The discretization errors are determined by the momenta 
carried by the quarks (and gluons) and are typically pro- 
portional to (a X momentum) 2 where a is the lattice spac- 
ing. The full mNRQCD action described in this paper 
has no tree-level O{a?k^)-eTY0YS, but has ©(a^a^fe^) er- 
rors due to radiative corrections. The same is true for 
highly improved light quark actions such as ASQTAD 
[IjiHalla] or HISQ [131 . Assuming that the constants of 
proportionality for the discretization errors arc the same, 
discretization errors are minimal if all quarks involved in 
the decay have momenta of the same size. 

The increase in discretization errors for the quarks in 
the B meson due to the boost of the momentum distribu- 
tion when going from zero velocity (7 = 1) to a non-zero 
velocity v is proportional to 

7'A|cD-A|cD- (1) 

Assuming that the quarks in the light meson share the 
momentum equally, each carrying momentum of order 
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FIG. 1: Estimate of optimal velocity, minimizing discretiza- 
tion errors, for B ^ F form factors (see text) as a function 
of for the tt, K and K* light mesons. 

Ppjl^ we expect that the increase in relative discretiza- 
tion errors for the light meson when going from zero mo- 
mentum to Pp is proportional to 

(^IPfI)'- (2) 

The total error is the sum of these terms with some co- 
efficients that we presume are of order unity. Noting 
that (j)B ^ Pf)'^ = (f' , we choose v to minimize the total 
error to give the optimal w as a function of . Inves- 
tigation with different reasonable choices of coefficients 
shows that the minimum error is compatible with setting 
the two above terms equal. The result is plotted in Fig- 
ured] for the TT, K and K* light mesons. Wc find that at 
maximum recoil a velocity of |d| « 0.7 would minimize 
discretization errors. Of course this is only a very crude 
estimate, and the optimal velocity depends on the details 
of the lattice computation. 

2. Statistical errors 

Lattice calculations are not only limited by discretiza- 
tion errors, but also by statistical errors. Unfortunately 
these increase exponentially when going to lower . Con- 
sider for instance the i3-meson two-point function with 
momentum p^, denoted as {B'' {pb,0)B{pb,t)) , which 
is required in the form factor computations alongside the 
pion two-point function and the B F three-point func- 
tion. The variance in the correlator is [28| 

a'ir) = {[B^{pB,0)BipB,r)][B^PB,0)BipB,T)]^) 
^\{BHpB,0)BipB,T))\\ (3) 

The correlator in the first line of ^ couples to the com- 
bination of a heavy-heavy (HH) meson at rest and a pion 
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at rest, so for large Euclidean time r, it will decay like 
exp(— (Mjjj/ + M^)t). However, the second line is simply 
the square of the two-point function, which will decay like 
exp(— 2i?B(p^)T) where Eb{Pb) is the energy of a B me- 
son with momentum pg. Since Mhh + -^tt < '^Eb{Pb)i 
the variance will be dominated by the first line at large 
r. This means that the signal-to-noise ratio approaches 
zero exponentially fast in Euclidean time r. 



— cx e ^ ^ 

cr(r) 



2 Mffff-iAf^^i 



(4) 



and at fixed r it decreases as the momentum pg in- 
creases. A similar analysis can be performed for the 
B ^ F three-point function and for the light meson two- 
point function. At lower q^, the momenta pg, pp and 
the corresponding energies are larger and hence the sig- 
nal decays faster, while the variance is independent jaf . 
(For an example with heavy-light correlators, see [29].) 

The above argument illustrates that using mNRQCD 
to extend the kinematic range of calculations requires the 
efficient use of techniques for reducing statistical noise. 
Already progress has been made reducing statistical er- 
rors using stochastic sources [sO]. Nevertheless, calcu- 
lations at lower q^ will undoubtedly require increased 
computational effort. In view of the opportunity for 
rare b s decays to discover or further constrain non- 
standard Model physics, via B K*-f, K^*H+l- for 
example, such effort is worthwhile. 



3. Heavy-quark expansion of the current 

Even in continuum mNRQCD systematic errors for 
heavy-to-light decays increase when going to lower q^, 
since the convergence of the heavy-quark expansion for 
the current mediating the decay gets worse. The heavy- 
quark expansion requires that all momentum scales for 
the light degrees of freedom are small compared to the 
mass of the heavy quark, which is approximately equal to 
the mass of the B meson, Mb ■ In the low-recoil regime, 
the only relevant scale is Aqcd ^ Afs, but for large re- 
coil the momentum of the light meson in the B rest frame 
is large. 

The light meson energy in the B rest frame can be 
written in a Lorentz-invariant way as 



E 



F.O 



PB-PF _ Ml + Ml - g2 



Mp 



2Mb 



(5) 



The light meson momentum in this frame is then \pp q | 



Ml. In Fig. O we plot the ratio \ppq\/Mb as 



a function of q^ for for the tt, K and K* light mesons. 
This ratio becomes almost 0.5 at q^ = 0, which has to be 
compared to Kqqd/Mb « 0.1 in the low-recoil limit. 



a. 0.2 




FIG. 2: The ratio \pp Mb as a function of for the n, K 
and K* light mesons. 



III. DERIVATION OF mNRQCD 

A. Continuum mNRQCD 

To derive the mNRQCD action in Minkowski space, we 
work in two frames, the optimal frame with coordinates 
X and the rest frame of the B meson with coordinates 
x' . The two frames are related by a Lorentz boost with 
velocity t>, 

X = Ax'. 

For the explicit form of A = A{v), see Appendix \X[ We 
denote the physical (full QCD Dirac spinor) heavy quark 
field in the two frames by \E'(a:) and '^'{x'). They are 
related by the spinorial representation of the boost, 

^(x) = S'(A)*'(a;'), 
'^(x) = W{x')S{A) . 

The spinorial boost matrix S{A) is defined in Appendix 
lAl The Dirac Lagrangian for is 



C'{x') = 4"(a;')(«7 • D' ~ m)^'{x'). 



(6) 



(The hat simply distinguishes a Dirac spin matrix from 
the 7 of the Lorentz transformation. Our convention 
for these matrices is given in Appendix [XI) Since the 
heavy quark is approximately at rest in this frame, we 
can approximate this Lagrangian very well by the stan- 
dard NRQCD Lagrangian. One approach to constructing 
this Lagrangian is by writing down all possible operators 
that are allowed by the symmetries of the theory This 
approach is described for example in [ij] and [3l| and 
has the advantage that it includes operators which only 
arise at higher loop order as, for example, four-quark op- 
erators. By matching to full QCD one finds, however, 
that these are suppressed by and will play no role in 
our analysis. 
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1. FWT transformation 

We use a Foldy-Wouthuysen-Tani (FWT) transforma- 
tion to derive the Lagrangian order by order in 1/m via 
field redefinitions, since this automatically generates the 
correct tree level coefficients of all operators. For a de- 
tailed description of the method, see (33 |. The transfor- 
mation can be written as 



in (jlOp . This is in contrast to the standard continuum 
"moving HQET" Lagrangian. 

Under the change of coordinates x = Ax' , derivative 
operators in the Lagrangian and FWT transformation 
transform like 

D'^^A\D,. (11) 
The transformation law for the gluon field strength ten- 



*'(x') = t;. 



(7) 



which defines the transformed field 4*'. (A corresponding 
transformation defines ^''(x')). The factor e^*™^ t' re- 
moves the additive mass term from the Lagrangian and 



^FWT is given by 



"^FWT exp 



X exp 



X exp 



1 
2m 
1 

1 



(ij ■ D') 



X [l + 0(l/m4)] . 



(8) 



(The chromoelectric and chromomagnetic components of 
the gluon field strength tensor are defined by — Fok, 
Bj = —\ejkiFki in Minkowski space). The resulting 
Lagrangian is 



D' g 

h — - 

2m 2m 



SB' 



9 



(^D""^-E' + ■ {D'xE'-E'xD'] 



(9) 



with 



F'^^ix') ^ A" ^A\F,^{x) 

leads to the following transformation for the chromoelec- 
tric and chromomagnetic components: 

E'{x') = -/(^E{x)+vx B{x) - ;^^'^('^ ' E{x))^ , 
B\x') = 7^B(a;) -vx E{x) - ::^v{v ■ B{x)) 



(12) 

Using (Uni), HT]) and the Lagrangian dH) can be 

expressed entirely in the new frame with coordinates x. 
Note that Lorentz invariance can be used to simplify the 
transformation in the following way: — u' ■ x' — u ■ x, 
where u' ~ (1,0) and u = (m°,m) = (7,71)). Similarly, 



D' 



D' 



D and 



{u ■ oy 



D^. The 



term with the adjoint derivative of the chromoelectric 
field can be written as D""^ E' = Dfu^F'"'. The other 
occurrences of the field strengths are simply replaced by 
([T^ . but we will not insert this expression explicitly for 
the sake of legibility. The Lagrangian becomes 

{u-Dy-D^ 

— I— 

2m 



£ = 



Z7°u ■ D + 



2m 



9 



8m2 



7" {Dfu^F^'' + ie,uiT.= A\ {D^, E[}) 



(13) 



cr^ 



Note that in ^ the adjoint derivative D''^'^ acts on 
E' only, whereas the standard derivatives D' act on all 
quantities to their right. 

As a result of the FWT transformation, all operators in 
the new Lagrangian commute with 7°, that is, the quark 
and antiquark components are decoupled to this order. 

The next step is to re-express the Lagrangian, ([5]), in 
terms of quantities in the frame x (which we will put onto 
the lattice). To this end, we define a new field "^[x) via 
the trivial transformation law 



^{X) EE ^'{x') 

¥(a;) = 'W'{x') 



(10) 



Note that in order to preserve the commutativity with 
7° we do not include the spinorial boost matrix S{A) 



2. Removing time derivatives in the Ham,iltonian 

Note that the operators of order 1/m and 1 /m? in (fT3|) 
now contain time derivatives. In the following, we will 
show how these can be removed via further field redef- 
initions to ensure that in the lattice computations the 
propagator can be obtained by solving an initial value 
problem using a time evolution equation. 

It is convenient to write the Lagrangian (|13p in the 
following form, 



£ = 7* 

with 

00 = 

01 - 
O2 



Oo + — Oi 

7m 



1 



(7m)' 



:02 



§ + 0(1/7713), (-;^4) 



2 



1 



§7 7° {Dfu^F'^'^ + ie,M^^A\ {D,, E[}) 
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We start by removing the time derivatives in Oi . To see 
how this can be done, we note that any field redefinition 
of the form 



exp U ^-(1), 

\7m ' 



^' = ^'(1) exp U 

\7TO 



will result in 



7TO 



On 



(1)1 



(7m) 



2^(1)2 



*(1 



(1) 



'C = 7*(i) 

with the new operators 

0(1)1 = Oi + {U, Oo}, 



0(1)2 = 02 + {C/, Oi} + C/OoC/+ 2 {'^"^ Oo} 



(15) 



Thus, we need to write Oi = 0(i)i — {U, Oq} with some 
operator U such that 0(i)i does not contain time deriva- 
tives. This is indeed possible: 



= On) 



1 



^ v ' 

= - {U,Oo} 

and we can now read off the operator U : 



U = [(7' - 1)^0 + (7' + 1)^^ • D] . (16) 

The next step is to remove the time derivatives (other 
than the adjoint time derivative, which acts on the gluon 
field strength only) in the new operator 0(i)2, given in 
psp . Similarly to before, we use a field redefinition 



*(i) 



exp 



(7771)2 



*(i) = *(2)exp 



(7777)' 



(17) 



now with an extra power of 1/(7777), so that the lower or- 
der terms are unaffected. The derivation of the operator 
V is given in Appendix IB] 



mNRQCD Lagrangian 



Finally, we rescale the fields 

*(2) = 



n/7 

*(2) = 



(18) 



to remove the factor of 7 in front of C. We arrive at 
the following result for the tree-level moving NRQCD 
Lagrangian in Minkowski space: 



if'Do + ifvD 



(v-D) 



27m 



9 



27777 



-JL^f (d""^ -E-v- {D'"'^ X B 



*^ fT,-{DxE'-E'x D) 

f{v-D, S • (v X E')} 



877772 



8(7+1)77^2 
(2-^^)5.0 



167772 
^9 .^0 



7" [Dl^'-v- D""^ ]{v-E 



4^2j^2 
3 



+0(1/777^). 



(19) 



As before, all terms commute with 7". We can therefore 
introduce 2-component fields ipvix) and ^i,(cc). 



Ipv 



to explicitly separate the Lagrangian into the quark and 
antiquark pieces: 



Hi 



iDo + iv-D - 

iDo + iv-D - 
+ 0(1/7772). 



(v-D) 



27777 

_ (t;.£))2 



27771 



27771 



2jm 



-a B' 



r B' 



(20) 



Terms with odd powers of 1/m (i.e. those without a fac- 
tor of 7*^ in ()19|) ) appear with the opposite sign in the 
antiquark Lagrangian. 

Note that we have chosen a particular notation con- 
vention for the 2-component antiquark field: ^t, creates 
an antiquark whereas ipv annihilates a quark. While the 
quark and antiquark terms in (|20p take a similar form, 
dictated by charge conjugation invariance, it should be 
borne in mind that ip^ and ^„ have these different inter- 
pretations when constructing the heavy quark and anti- 
quark Green functions. As an aside, note that our new 
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QQ rest frame lattice frame 



D' 

gE' 

gB' 


2 

mVral 
2 4 

m v^^i 


Dt 
D 
gE 

gB 


7m(llj(,l + VVrcl) 
7m(Wrcl + UUrol) 


77Tl?;rcl 
7mWrcl 


D'D' 
u'-D' 


2 2 


DD 
u-D 


2 2 / / 2 2 2 \ 
"1 ^^rcl {^1 m «rcl) 


2 2 



TABLE I: Power counting rules appropriate for mNRQCD 
with heavy-heavy mesons. In the large velocity limit (last 
column), the Lorentz boost removes the differences in order 
found for NRQCD, giving A ~ -D and ~ B. In the last 
two rows note that the naive power counting rules can give 
the wrong counting (see text). 



result differs sli ghtl y at order \/m? from the one 
given in Refs. H,!!!!!!. 

Let us now summarize the tree-level relation between 
the full QCD field *(x) and the moving NRQCD two- 
component fields V'f (2;), £,v{x)- 

^{x) = S{K) T,wT e— A^^ ^ ^ (21) 

where TpwT is the FWT transformation ([8]) expressed in 
the frame x, i.e. 

Tfwt = exp ( I exp ( '^T'^.T I X ... 



2ra 



\ (2m)2 



and 



removes the unwanted time derivatives in the Lagrangian 
{U and V were defined in equations (fT6|) and (jBSp . re- 
spectively). 

The field redefinition (pij) can be used to obtain tree- 
level expressions for currents containing the heavy quark 
in calculations of decay constants and form factors, as 
discussed briefly in section IVD] 



B. Power counting 

When deriving the mNRQCD Lagrangian in the pre- 
vious section, we were formally expanding in powers of 
1/m. As is well known from heavy-quark effective theory, 
for heavy-light systems such as B mesons, the expan- 
sion really is in Aqcd/"i with the QCD scale Aqcd ^ 
500 MeV. The Lorentz transformation docs not affect the 
power counting and thus the Lagrangian (|19|1 is complete 
through order (Aqcd /to)^. For heavy- heavy mesons such 
as the T, the situation is more complicated. In the frame 



where the meson is at rest, the power counting is gov- 
erned by powers of v-^cii the small non-relativistic inter- 
nal velocity of the heavy quarks inside the meson p^ . 
For T systems, one has v'^^-^ ~ 0.1. It turns out that all 
terms of the Lagrangian ([9]) are of order v^^y or lower, 
but one term of order v^^y is missing. By expanding the 
expression for the relativistic kinetic energy in powers of 
the residual momentum fc, 



k" 



2m 8m^ 16m^ 

and replacing fc by the operator iD, we see that we 
must include the operator D'^/iSnr^) into ^ in order to 
obtain accuracy to order vf^^. The corresponding term in 
the moving NRQCD Lagrangian can be obtained in the 
same way. 



Ekin — I'm = \J {■ymv + fc)^ + — 7m 
1 '.2 



— V - k 



1 



1 



(k'-iv-kr) 

-{v-k, k^} + 2{v-kf) 

-fc* + 3{fc^ {v-kf} 
-5iv-k)^) 



+ ... 



Thus, the operator 
1 



(r>4 - 3 {D^ {v -Df}+ 5iv - Df) (22) 



must be included into the moving NRQCD Lagrangian 
(fT9)l . We ordered the terms with products of {v - 
D) and in the form of anticommutators, as the 
anticommutator-ordering is what one would have ob- 
tained from field redefinitions. 

For heavy-heavy mesons at 1; = 0, the power count- 
ing is different for temporal- and spatial components of 
Lorentz vectors but they will mix in a frame with v ^ 0. 
The rules for both v = and v ^ are summarized in 
Table H 

Care has to be taken when dealing with quantities like 
D ■ D and u-D; their power counting cannot be derived 
by naively multiplying the power counting rules for each 
factor. For example, for v 1 the product D - D does 
not scale like {'jmv^ci)'^ but as m^v^^^ instead. The correct 
values are shown in the last two rows of Table HI 



C. Euclidean mNRQCD 

The Euclidean action Se = J (i^XE .C-EixE) can be ob- 
tained from the Minkowski space action S = J d'^x C{x) 
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in the usual way by making the formal replacements 



A{x) 
Aq{x) 



A{xe)-, 
iAi{xE), 



A — 



'IT, 



so that the integration measure and derivatives become 
d^x — !■ {—i)d'^XE, do ~^ idi. Finally, the result must be 
multiplied by (— i). In the following, we drop the sub- 
script E ("Euclidean"). Note that we do not introduce 
Euclidean gamma matrices in this paper; the same defi- 
nition as in Minkowski space is used (see Appendix |^ . 

It is also convenient to define the relation between the 
chromoelectric field E and the 4-dimensional Ff^i, with 
a different sign in Euclidean space, i.e. Ej = Fj4, while 
the definition of the chromomagnetic field is unchanged, 
Bj = —^ejkiFki. 

With this definition, turns into the symmetric 

form 



E' ^-/(E + iv X B 



7 



j+1 



-v{v ■ E] 



B' = -/[B + ivx E- 



7+1 



-v(v ■ B 



(23) 



The Euclidean Lagrangian, in which we now include the 
relativistic correction term (|22p . becomes 



^"Di - i^^v-D - 



2jm 



9 



2jm 



SB' 



7° {{v-D, D^] - 2{v-D) 



8^2' {^D^'■E + v■{D^'xB) 

— ^7° ^ ■ {D X E' - E' X D) 

87771^ 

^ ^f{v-D, ^-{vxE')} 



^7° 



8(7 + 1)to 

ic 2 ^ 1^4 



%v ■ D""^ ](v-E) 



4^2^2 
1 



f{v-D, T,-B') 



+b{v ■ D) 



(24) 



As in (PO)) one can introduce two-component fields for 
quark and antiquark. It turns out that in Euclidean 
space, the antiquark action can be obtained from the 

quark action by replacing — > (Ci) , V4 ~* {£.v)^ , 
V {—v) and taking the complex conjugate of the whole 



action kernel. This is an important result, because it im- 
plies that the Euclidean antiquark Green function can be 
obtained from the Euclidean quark Green function in a 
frame with the opposite boost velocity, —v. We define 

g''^'"\x, x') = {iv{x)^l{x')). Writing out color, spin 
and position indices explicitly, one then has 



G 



(x, x') 



G 



G 



i-v) 



c s cs 

t 



(a;', x) 

S 

(x', x). 



(25) 



IV. LATTICE mNRQCD 

A. Construction of the Hamiltonian 

We construct the lattice moving NRQCD action such 
that for u = it reduces to the previously used lattice 
NRQCD action with conventions as in [3J|. Thus, the 
quark action has the form 

S^^=Y,^li^^^)[Mx,T)-K{T)^,{x,T~l)\ (26) 
with the kernel 



K{t) = [I- 



SH\r 

2 

^^olr-l 

2n 



(27) 



Note that the heavy-quark Green function for the action 
([26|) satisfies the evolution equation 



{x, T, x', t') = K{t) {x, r - 1, x' , t') 



(28) 



For this, it is crucial that the Hamiltonian does not con- 
tain time derivatives (other than the adjoint time deriva- 
tive of the chromoelectric field). 

This split into leading-order kinetic terms Hq and 
higher-order corrections 5H which satisfies time-reversal 
symmetry was introduced in [l4| . Other than consistency 
with previous work, there are no strong arguments (such 
as computational load, numerical stability or size of dis- 
cretization errors) for the relative ordering of Hq and &H 
in the action. The time derivative in Ij26p is implemented 
as a backward (rather than forward) difference operator 
as this prevents mean-field corrections to the wavefunc- 
tion renormalization 14]. 

The leading evolution due to from one lattice time 
slice to the next is effectively divided into 2n smaller 
steps to avoid the well-known instability in the discretiza- 
tion of parabolic differential equations (see, for instance. 
Sec. 19.2 of Ref. [1^). In this way, one can allow the 
highest momentum modes in the theory to come into 
equilibrium, while avoiding the need for a very small lat- 
tice spacing which would render the theory too expensive 
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to simulate. For NRQCD, where Hq is always positive, 
the integer-valued stability parameter n has to be chosen 
such that 



1 - 



Ho 
2n 



< 1. 



(29) 



In the free field case this condition can be satisfied by 
choosing n > 5/ {2am), and gluons are known to reduce 
the factor of 3/2 shghtly [H. 

In moving NRQCD, Hq can be negative for values of 
k pointing opposite to the frame velocity. In this case 
the two-point function will grow exponentially, but this 
is physical as we find the same behavior in the contin- 
uum. In our numerical simulations, which included boost 
velocities up to w = 0.6, we did not encounter any insta- 
bilities with n = 2, am = 2.8. 

The lattice Hq and SH are defined as 



Ho = 
6H = 



A(2) 



^(2) 



2jm 
i 



2jm 

a B' 



(30) 



3) 



47-^™ 

g 

'8m2 
9 



(^i{A^ ■ E - E ■ A^) 



V ■ (A^^ X B 



A"^ X E' 



E' X A' 



9 



8(7-1- l)m2 
{2-v^)g ( 



{« 



. ad 



16to2 
«5 



A^, o- • (t> X i;')| 



, ad 



B 



1 



-^((a«)'-.{a<^a?.} 



5Al^) 



(31) 



The lattice derivative operators and field strength arc 
defined in Appendix [C] Note that in the continuum the 
Leibniz rule D^'^-E = D E~E D holds. For consistency 
with previous work we discretize the right hand side of 
this expression on the lattice. However, the other adjoint 
derivatives in the action, which enter only at v 7^ 0, are 
discretized as lattice adjoint derivatives. This is more 
efficient and for the term {v ■ E) it is crucial since it 
avoids a time derivative acting on the quark field. 

Note that in the static limit (m 00) one has Ho = 
—iv ■ A^. The symmetric derivative A^ leads to zero- 
energy modes at the corners of the Brillouin zone ( "dou- 
blers"). With a finite mass, these doublers are shifted to 
higher energy due to the second-order derivatives in Ho- 
However, the second-order derivatives are suppressed by 
a factor of 1/ (27m) and hence 7m must not be too large. 

The terms in SHcort provide the spatial and tempo- 
ral lattice spacing improvement. We perform tree-level 



Symanzik improvement to order 0(0^), as explained in 
the next section. This means that the we expect the 
leading errors to be of order 0{asa'^). 



1. Improvement corrections 



An C(a'*)-improved version of Ho is given by 



Ho 



= —iv ■ A — 



^(2) 



(32) 



with the improved derivatives given in Appendix [Cl 
However, we do not simply replace Ho by Ho- Let us 
first consider the time derivative in the lattice action. 
Improving it in the standard way would introduce next- 
to-nearest neighbor couplings, preventing the use of an 
evolution equation like (|28p . Instead, we try to find an 
operator Hq such that (explicitly re-introducing the lat- 
tice spacing a) 



1 - 



oHl 
2n 



exp 



(-i 



Hr 



(33) 



which yields a more continuum-like behavior ^\^. We 
obtain 



aHX 



2n 



1 — exp 



aHp 
2n 



(34) 



One could now replace Hq Hq in the lattice action. 
However, for performance reasons and consistency with 
previous work, we choose to put all correction terms into 
6H. We consider the operator on the right-hand side of 
the temporal link in the lattice action (|27p : the operator 
acting in the timeslice at time r — 1. Then SHcon, the 
lattice spacing improvement term in (|3ip is defined by 



1 - 



a_Hl 
2n 



= 1- 



oHq 
2n 



1 - 



(35) 



for SHcoii- This gives 



a SHc. 



= 2 



= 2 



1-1- 



1-1- 



aHp 
2n 

aHp 
2n 



exp 



oHl 
2n 



aHo 
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and, expanding in powers of a, 



Hi - nHl + 2nHoHo 



24n2 



- {2 + 3n + n^)H^ 

+ {3n + 3n^)HlHo 
-Sn^HoHl + n^H^ 



192^3 



+ Oia^). 

The term C = Hq — Hq is of third order, while Hq is of 
first order. Neglecting all operators of order 5 and higher, 
we obtain 



a 6Hc 



aC- — {Hl + n[C, Ho] 
An 

aJ^Hl (2 + n)a^H^ 



12n2 



(36) 



Had we considered the operators on the left-hand side 
of the temporal link in the lattice action (|27p instead, 
the ordering of Hq and Hq would be interchanged, and 
this would change the sign of the commutator [C, Hq\ in 
pop , thereby cancelling the term in the lattice action up 
to operators of order 5 and higher. We therefore remove 
this term on both sides. 



Let us go back to lattice units now. Writing Hq 
A + B with 

A = -iv ■ A^, 

A(2) - Ai') 



B = 



27m 



we obtain 

^Hrnrr 



Ho -Ho- — {A^ + {A, B} + B^) 
An 

{A^ + {A^ , B}+ABA) 



Ylv? 



.(2 + n)^4 
64n3 



(37) 



For performance reasons, we replace some 3rd- and 4th- 
order derivatives in p7p by more local expressions (the 
resulting change is of order 5 or higher): 

(t;.A±)3 ^ A(3), 

{^.A±, A(2)} ^ 2A(3), 

(t;.A±)4 ^ aW, 

(A(2))2 _> 

{(^.A±)^ A(2)} ^ {A(2', A(2)}, 

{(^.A±)^ A?)} 2AW, 
(t;.A±)A(2)(^.A±) ^ AW, 
(v A±)A(2)(^;. A±) ^ i(7;- A-)A(2)(t;. A+) 

i(^;- A+)A(2)(i;. A"). 



This finally gives 



5Hcorr — Ho — Ho 

If ^ {z^.A±, A(2)}-2*a13) (A(2))2-{A(2), A?)} + a1^) 



in \ 2jm A'^'^m? 

1 ^ ^ {A(2), Ai')}-3Ai^) + i((t;.A-)A(2)(^.A+) + (^.A+)A(2)(^.A-)) 



12n2 V " ' 27TO 
(2 -I- n) 



64n3 



AW. (38) 



The result (j38p can be simplified further since most op- 2. Radiative corrections 

erators are already in the Hamiltonian. 

In principle, all operators in the Hamiltonian are mul- 
tiplied by coefficients Ci which contain radiative cor- 
rections that correct for lattice artifacts appearing be- 
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yond tree-level, including the missing contributions of 
UV modes with momenta greater than the lattice cut- 
off; \k^\ > n/a. They can be expanded as a power series 



(0) 



+ ascf^ + . . . + (a.)"cf ) + 



where the tree level c, 
„(") 



(0) _ 



1 and the radiative correc- 



tions c\ depend on the bare quark mass and the frame 
velocity. These radiative corrections are calculated us- 
ing lattice perturbation theory by matching standard on- 
shell processes computed in mNRQCD with the contin- 
uum counterpart. Four-quark operators can only arise at 
0(0^) and for this reason will not be considered in our 
analysis. 

For the calculations in this paper, we use the tree level 
values of the couplings a. However, we account for a 
large amount of the expected renormalizations via tad- 
pole improvement. 



B. Tadpole improvement of the Hamiltonian 



The tadpole improvement of perturbation theory is dis- 
cussed in subsection IV HI 



V. RENORMALIZATION OF mNRQCD 

In the previous sections we derived the tree-level con- 
tinuum mNRQCD Lagrangian and its lattice version. 
The radiative corrections to the couplings Ci include a 
renormalization of the external momentum whose origins 
are discussed below. The momentum renormalization is 
important because it is the coupling of the v ■ D term (= 
Pq ■ k/jm, Pq = jmv) in the action which is leading or- 
der in the l/m expansion. The momentum renormaliza- 
tion must be well-determined for accurate results. Fortu- 
nately, as described below, approximate reparametriza- 
tion invariance ensures that this renormalization is small; 
the renormalization constant is close to unity. 



Derivation of the mNRQCD renormalization 
parameters 



It is well-known that the perturbative expansion in the 
bare lattice coupling is poorly behaved. Tadpole dia- 
grams, which do not contribute in continuum schemes, 
give large contributions to coefficients multiplying powers 
of the bare coupling. Tadpole improvement (also known 
as mean-field improvement) fixes this problem by resum- 
ming diagrams containing tadpoles [36|. As tadpole im- 
provement reduces the size of perturbative corrections, 
even the tree-level couplings in the action will give accu- 
rate results. Gauge links Ufj, and U"^ in the action and 
operators are divided by a factor uq which is designed 
to correct for the fact that the expectation value of the 
mean link (using some gauge-fixed or gauge-independent 
definition) is much less than unity. We choose uq to be 
the mean link in Landau gauge. The fourth root of the 
mean plaquette is another frequently used definition of 

Uq. 

Care has to be taken when replacing 1-^ U^/uq 
UJ^/uq in the action. The action is composed 
of Wilson lines or "paths" . If, due to application of a 
lattice derivative for example, the product U ^{x)Ujj^{x) 
appears, one should not multiply by a factor of 1/uq 
since the product is trivial and does not contribute to 
tadpole contamination. Some paths are not explicit in 
our simulation code, where we evolve the heavy quark 
green function by subsequently applying the individual 
blocks of the action kernel (P7)) rather than expanding it 
in terms of paths first. Explicitly coding (|27p in terms 
of products of link variables would be forbiddingly time- 
consuming. Therefore we only take into account link-pair 
cancellation separately within Ho and 5H . Also, no ex- 
tra cancellations are made when derivative operators act 
on field strengths in 5H . 

For perturbative studies the tadpole counter-term 
must be computed to the appropriate order in 0{as). 



and Ul 



The low-momentum properties of the moving heavy 
quark inverse propagator can be expressed as a general 
power series in the energy p4 and the three-momentum 
p. The coefficients of this power series determine the 
renormalization of the wavefunction Z^, the quark mass 
Zm, the shift in the origin of energy Eq and of the frame 
velocity Zy. 



1. Wavefunction renormalization 

The wavefunction renormalization Z^ can be com- 
puted using the following simple arguments. The tree- 
level quark propagator is given by: 



Go(z) 



Z- 



where z = e*''"' and 



zo = 1 - 



go(p) 
2n 



1 



(39) 



(40) 



Then z — zq \s the on-shell (tree-level) value. At one 
loop 



G"i(z) = Go\z)-a,S(z) = Z^ 



_lZ- Zl 



where as^{z) is the self-energy (to order as), containing 
both rainbow and tadpole diagrams. Let the new "one- 
loop" on-shell value be zi , which is the solution of 



G-i(zi) = Goi(zi)-a,I](zi) = 0. 



(41) 



Expanding Y,{z) around the new on-shell value we have: 



E(z) = S(zi) + (z-zi) |- 



(42) 
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Therefore 



1 



G ^{z) = -(z-zo)-as S(zi) 
z L 

+ (z - Zl) — 

oz 



Eliminating zq in this expression in favor of zi using (|4ip , 
we obtain 

G-\z) = i(z - zi) ( 1 - a, [e(zi) + z ^ 

(43) 

Thus, as zi — zq CI {as), the wavefunction renormahza- 
tion is, at one loop, 



Z^h = 1 + as 



1 + as 



S(zo) + z — 
oz 



as 

dp4 



on shell 



and the renormalization constants can be expressed in 
terms of the coefficients S^- in the expansion 



We find 

Ea - 

Zm - 



(0) 

l-a.(EW+s(«)), 

l + a.((s(^)+I]l"V7V(E(°)+E«)) 



(48) 



and have for the renormalization of the external momen- 
tum P = jruirVr = ZpPo with 



Zri 



l + a.(I](")-E(°)). 



(49) 



^44^ In actual calculations we consider the real parts of pa- 
rameters Ti^^\ It is convenient to define 



2. Other renormalization parameters 

To derive the other renormalization parameters, we use 
the following argument which can easily be extended to 
higher order kinetic terms [37| . At tree level we have in 
momentum space (up to 0{p^)): 



— {v ■ pY 
27771 



Ho{p) = V p 
8H(p) = -l-{vpf + 



(45) 



By combining this with (|40p and expanding in p we find 
that the pole in the tree level propagator ([55]) is given by 



UJ = ujo{p) = vp 



p^ — [v ■ pY 
27777 



(46) 



where uj = —ipi is the energy in Minkowski space. At 
one loop the inverse propagator is 



G(p,c.)-i = l-e--o(P)-a,I](p,c.o(p)) 



so that 



uj{p) = tJo(p) - asS(p,cjo(p)) 

- (VR ■ pY 



(47) 



vr-P 



2^RmR 



- as6uj{p) 



with Vr = ZyV, ^R = {1 — v^)"^/^, rriR = Zmm and 
asSuj{p) — Eq + . . . . Here and in the following we as- 
sume that the boost velocity points in one of the lattice 
directions, which guarantees that only the magnitude of 
V is renormalized. The self energy can now be expanded 
in small momenta 



I](p, u) = So(w) + J:yiuj)v-p + Ei(w) 



P 



2'ym 



no 
ill 

^2 



Re = E(0) 



-ReE, 



(1) 



Im 



9E 



p=0 



Re e'"^ = 777iRe ^ 
dpi 

ReE(«) = iRe |^ 
V dpx 



92 E 



p=0 



(50) 



p=0 



taking the frame velocity v to lie in the x-direction. The 
renormalization parameters are then expressed as 

Z^ ^1 + as{Vlo + rii), 
Zy = 1- as{Vly - fii) . 

z,„ = i + as[n2- ni) + asiVLy - ni)v'^-f'^ , 

Zp = I - asiSly - ^2) . (51) 

B. Dispersion relation and energy shift 

The renormalized dispersion relation in (I47p has to be 
compared to the corresponding expression in QCD 



^(QCD)(p) ^ y'(^^^^^^+p)2+^2^ (52) 

= iRniR + VR-P 

^ P^ - {vnrP)^ ^ 
'^iRrriR 

from which one obtains a shift in the zero point energy 
of a heavy quark of 

Cy = w(Q^°)(p = 0)~L^(p = 0) (53) 

= jRrUR + Eq . 
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We write = 7m(l + as5Cy -\- 
correction is given by 

5Cv ~ — f^i 



. ) and the one-loop 



7m 



(54) 



The shift and the renormahzation of the external 
momentum can be obtained nonperturbatively by com- 
puting the energy (p) of a heavy-heavy system which 
is up to lattice artifacts given by 



ML- 



A corresponding dispersion relation with 2Cy 1— > C„ and 
2Zp i—f Zp holds for heavy-light mesons containing only 
one heavy quark. We will compare values for the energy 
shift and the renormahzation of the external momentum 
calculated in perturbation theory and nonperturbatively 
using the dispersion relation of heavy-heavy and heavy- 
light mesons in section IVlIl 



C. Reparametrization invariance 

One thing we expect from our results is that, because 
of lattice reparametrization invariance ,33] , the deviation 
of the momentum renormahzation parameter Zp from its 
tree level result is much smaller than for other renor- 
mahzation parameters. Reparametrization invariance is 
a symmetry that has been studied in the context of 
heavy quark effective theories [sl, [H, H^]. This sym- 
metry arises from the fact that the division of the full 
momentum p into a "fixed" external part m u and a "dy- 
namic" residual part k is not unique. We can always 
write p = m u + k ~ m u' + k' where k' ^ k — me, 
u' = u + e. The 4-velocities u and u' have unit norm 
which implies the constraint on e that 2e ■ u + e ■ e ^ 0. 
It can be shown [38l . [3^ that this reparametrization of 
the full momentum is a symmetry of the effective heavy 
quark Lagrangian in the continuum. 

Because mNRQCD is a non-relativistic formulation 
Lorentz symmetry is not manifest in the action. This 
is apparent from: the form of the FWT transformation 
([5]); the field redefinition f section [ill A 2p required to re- 
move time derivatives in the Hamiltonian; the trunca- 
tion of the action to a given order in 1/m; and the non- 
relativistic field normalization p^ . To adapt the dis- 
cussion of reparametrization invariance to mNRQCD we 
study the ambiguity in division of the total 3-momentum, 
p = ^{v)m V + k, keeping fixed since it is a parame- 
ter in the Hamiltonian. We first consider a simple action 
with Hamiltonian 



Hn = -iv ■ D - 



27m 



(55) 



omitting the term {v ■ £))^/(27m) for the moment. This 
action is invariant under the transformation 



tjme-x 



(56) 



with 2v ■ e + e ■ e = 0. This constraint ensures that 
\v'\ = \v\. This is an exact symmetry which implies that 
the external momentum Pq — jmv is not renormalized 
as the relative coefficients of the two terms in ([55]) are 
fixed even after renormahzation. 

On the lattice, where we use the discretized Hamilto- 
nian 



H, 



(lat) 



.± 



A(2) 

27am' 



(57) 



this symmetry is broken. Under (|56p Hq^*''^ transforms 
according to 



H, 



If V is chosen along a lattice axis Vj — v6ji, say, then us- 
ing the constraint on e the factor vjCj can be replaced 
by — i|ep(5ji which is small for small e. We might 
therefore expect the breaking of reparametrization in- 
variance by lattice artifacts in this case to be small. In 
the corresponding derivation with improved derivatives 
in Hq^^\ we find that the lattice artifacts which break 
reparametrization invariance are of 0{a'^). 

Reparametrization invariance is broken even for the 
continuum theory unless the FWT transformation and 
the truncation of the action as a series in 1/m respect it. 
The field redefinition, designed to remove time deriva- 
tives in the Hamiltonian, must also be invariant under 
the reparametrization transformation. This will be sat- 
isfied only if the velocity and the covariant derivative 
appear in the combination (ssl . [39| 



iD 

27m 



(58) 



This implies that terms of different order in 1/m are 
mixed by the reparametrization transformation and so 
any truncation of the action as a series in 1/m will 
break this invariance. It would be possible to include 
selected higher-order terms in 1/m by rewriting the ac- 
tion in terms of the combination (j58p but in practice 
this is unnecessary since the approximate reparametriza- 
tion invariance of the action is sufficient to restrict the 
renormahzation Zp of the total quark momentum Pq to 
be close to unity. It would also introduce extra terms 
of little significance in the non-relativistic expansion but 
which are expensive to evaluate computationally for the 
lattice theory. In any case discretization breaks the in- 
variance as already discussed. We shall compute Zp both 
perturbatively and nonperturbatively. 

The mixing is evident in our simple example above. 
It is easy to see that adding the term {v ■ D)"^ / {2jm) 
will break the invariance for non-zero frame velocities 
even in the continuum. This breaking is proportional 
to v'^/{2jm), so it increases to reach a maximum at 
V = -^72/3 « 0.8 and then drops to zero due to the sup- 
pression by 1/7. Numerically we find this behavior in 
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our perturbative results for the simple action we discuss 
in Appendix IE II The one-loop contribution to the ex- 
ternal momentum renormalization vanishes for small v, 
rises to a maximum at z; « 0.75 and then drops again. At 
this velocity we also computed 6Zp with the action (j55[) 
both with naive and improved derivatives. We find that 
the use of improved derivatives reduces dZp by roughly 
a factor 2. 

Our numerical results (see Table IIVI to be discussed 
in Sec. IVI|1 do indeed show that on the lattice Zp is very 
close to 1 for small frame velocities. For larger frame 
velocities the perturbative results show a deviation of Zp 
from the tree level value of at most 10% for practical 
choices of frame- velocity v. 



D. Current construction 



between on-shell states in the continuum and lattice the- 
ories. 

Note that the renormalization of the boost velocity also 
affects the spinorial boost matrix 5'(A). We have for bare 
quantities 



S{Aiv)) 



1 



I 1+7 J (T ■ V 

^2(1+ 7) \1(T V 1+7 



and the renormalized matrix is obtained from this by an 
additional Lorentz boost, 

SiA{vn)) - S{AiSv))SiA{v)). 

(No Wigner rotation is needed here as only the magni- 
tude of V is renormalized for v pointing in one of the 
lattice directions.) We find 



For calculations of hadronic matrix elements of weak 
interaction operators involving the heavy quark, the con- 
tinuum QCD currents must be replaced by appropriate 
lattice currents. Let us, for example, consider the vector 
current 

where q{x) is the Dirac field of the light quark and ^>{x) 
is the Dirac field of the heavy quark. At tree-level, it 
suffices to express '^'{x) via the Euclidean version of the 
field redefinition (PT|) . 

Recall that Eq. ^ contains a factor of e"™^'"'^" which 
removes the mass term from the Lagrangian. For a heavy 
quark, the lower two components of the non-relativistic 
field ^'{x') in Eq. ^ are zero, so that ^°^'{x') = ^''(a;') 

and hence e-™"^"'T'''^'(a;') = e-™"="'^'(x'). Since the 
FWT transformation in this frame does not contain time 
derivatives, the factor g"*™^ can be moved to the left 
of Tp^T. (In the antiquark case, where the upper two 
components of ^'{x') are zero, one has e+""^ .) 

Performing the other steps of the derivation in Section 
lllll again. it then follows that also the factor of e~"" "'^ ^ 
in Eq. (|2ip can be moved to the left of Tfwt in the case 
where £,vix) — 0. Thus, in correlation functions the factor 
g-im u-x 7 ^j-ivially shifts energy and momentum and can 
be removed. We obtain 



S{AiSv)) 



1 



^6v ■ <j 1 




J^^{x)^q{x)r^S{A)T,, 



For on-shell quantities, time derivatives in TpwT and A^^ 
can be eliminated using the equations of motion, Z?4 — 
iv ■ D + 0{l/m). The continuum derivatives are then 
replaced by lattice derivatives. 

Beyond tree-level, additional lattice operators are re- 
quired and matching coefficients must be introduced to 
correct for the different ultraviolet behavior of QCD and 
lattice mNRQCD. These matching coefficients can be 
computed perturbatively by comparing matrix elements 



with 



ov = = agOZyj V. 

1 — V ■ Vr 

We will not consider the current matching any further 
here; this will be discussed in another paper. 



E. Lattice Perturbation Theory 

Feynman rules for lattice actions are complicated and 
for all but the simplest cases an automated procedure is 
needed to obtain them. The formalism for this is due 
to Liischer and Wcisz [4l|. This was extended by Nobes 
and Trottier ^42,] and Hart et al. [H, H [H to include 
both relativistic and non-relativistic fermion actions such 



as HISQ [27[ and, as used here, mNRQCD. In this paper 
we use the implementation of Hart et al. [H, 13, |45i| to 
compute the one-loop self-energy I](z) for various choices 
of mNRQCD Hamiltonian. The Feynman rules, vertices 
and propagators, are generated in machine-readable form 
using the Python program HiPPy and then used in the 
Fortran 95 code HPSRC to construct the diagrams and 
carry out the loop momentum integrations. The latter 
are done using VEGAS [H, \^ or, in the case of small 
lattices, by mode-summation. All perturbative results 
presented in this paper are obtained on an infinite lattice. 

As more correction terms are added to the action, the 
number of terms in the perturbative expansion grows 
very fast, and so does computation time. We have used 
a version of VEGAS that has been adapted to parallel 
computing using MPI (Message Passing Interface). 

The diagrams we evaluate to obtain the heavy quark 
self energy at one loop are shown in Fig. [S] The renor- 
malization parameters require derivatives of the self en- 
ergy. The derivatives of the Feynman rules were calcu- 
lated exactly (rather than from small finite differences 
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FIG. 3: Diagrams to be evaluated: (a) rainbow diagram and 
(b) tadpole diagram. Numerical calculations show that con- 
tributions to the heavy quark self-energy from diagram (6) 
are approximately an order of magnitude bigger than those 
of diagram (a), demonstrating the crucial importance of tad- 
pole improvements for any lattice-based perturbation theory 
calculation. 



due to their associated errors and instabilities) and then 
automatically combined to form diagram derivatives us- 
ing code based on the TaylUR package [H, \^ (which 
overloads arithmetic operations so as to respect Leibniz's 
rule and the chain rule). 

As an alternative to perturbation theory based on loop 
integrals, renormalized quantities may be measured by 
simulation in the weak coupling regime of the theory (i.e. 
at high j3) [53,[5l| on small lattices using 't Hooft twisted 
boundary conditions [4l|, [s^l . While not the subject of 
this paper, knowledge from analytic calculation of the 
one-loop corrections allows accurate fitting to extract the 
two- loop contributions. To implement twisted boundary 
conditions is straight-forward; it requires the spectrum of 
the momenta used to be appropriately modified and the 
vertices to carry a momentum-dependent phase rather 
than the usual color factor. We will discuss such cal- 
culations for mNRQCD in more detail in a forthcoming 
publication. 



Contour shift 



For a Euclidean lattice field theory the energy inte- 
gral is nominally over the unit circle \{z = e*'^*)! = 1. 
However, the positions of the poles in the integrand are 
functions of the loop three-momentum and care must 
be taken that no pole crosses the contour: the contour 
must be distorted to avoid this happening. In partic- 
ular, the heavy quark pole must remain inside the 
contour of integration in order to represent a forward- 
propagating heavy quark. This can be done by choosing 
\z\ — R,R > 1 where R is chosen so that the contour is 
large enough to enclose Zh and as distant from any pole 
as is possible to improve convergence of the integration. 
In Fig.|4]we show the position of the poles in the z plane. 

This contour shift applies to the case of the rainbow di- 
agram Fig. [3^ but is not necessary for the tadpole graph 
in Fig. as the poles in the gluon propagator corre- 
sponding to solutions moving forward/backward in time 
always come in pairs with = 1. 
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FIG. 4: Position of poles in the complex z plane and inte- 
gration contour (dashed circle). The two poles in the Wilson 
gluon action are z± with z+z^ — 1 whereas the heavy quark 
pole can be found at z^. If Zh > z- we shift the contour 
according to 2 i~+ Rz with R = ^ZhZ+ > 1. 



Finding the pole of the heavy quark propagator is 
straightforward as the Lagrangian only contains first or- 
der time derivatives . Exact expressions for the posi- 
tion of the poles of the Wilson gluon action can also be 
derived. These and the extension to more complicated 
gauge actions are discussed in Appendix |f1 There we 



show that Izll^P^I < z_ < 1 < z+ < |z!^'"p^| so that the 
contour shift derived for the Wilson action remains valid. 

The additional contour shift which is necessary when 
formulating the theory in Euclidean space has been dis- 
cussed in the literature [13, In Ref. [s^l, Aglietti 
et al. conclude that deriving Feynman rules for HQET 
in the Euclidean theory is problematic as a simple Wick 
rotation will generate unphysical solutions propagating 
backwards in time. However, in a subsequent paper [55| . 
Aglietti extends the analysis and realizes that this is 
due to an incorrect rotation of the integration contour 
to Euclidean time. To avoid crossing the heavy quark 
pole at i; • fe it is necessary to rotate the contour around 
—A = v-k — S instead of the origin of the ko plane (see 
Fig.ED. 



^(imp)| 



G. Treatment of infrared divergences 

To deal with infrared divergences, we note that any 
lattice theory has the same infrared behavior as the 
equivalent continuum theory. Therefore we consider the 
diagrams of Fig. [3] where lattice Feynman rules have 
been replaced by equivalent continuum ones (noting 
that the two-gluon vertex is still present in continuum 
(m)NRQCD). 

To analyze the infrared behavior of these diagrams we 
first perform the integration over the temporal compo- 
nent of the loop momentum as a contour integration, then 
look at the behavior of the remaining three-dimensional 
spatial integral for small loop momentum. 
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FIG. 5: Wick rotation to Euclidean space for continuum 
HQET; the integration contour is shown as a dashed line. If 
the heavy quark pole at ojh = v k lies to the left of the imagi- 
nary axis the contour has to be rotated around — A = v-k — 5. 
The gluon poles are denoted by ll}±. 



In the non-moving case {v — 0) this can conveniently 
be done in spherical polar coordinates; for the moving 
case we need to take into account the fact that the ex- 
ternal velocity introduces a preferred direction. It is con- 
venient to take the velocity to lie along the x-axis, for 
instance. 

After performing these calculations we see that the 
rainbow diagram Fig. [3^, as well as the tadpole diagram 
Fig. [3Jd and all derivatives of the tadpole diagram are 
infrared-finite; however the derivatives of the rainbow di- 
agram behave for low momentum as ~ / ^ and thus are 
logarithmically divergent. To regulate this divergence we 
introduce a small gluon mass A, which we may do because 
the rainbow diagram has Abelian color structure. 

To find the infrared behavior of fii , and Qy we per- 
form the analytic calculations as detailed above, keeping 
track of all prefactors in the integration. After doing this, 
we obtain the infrared-divergent part of the derivative of 
the rainbow diagram: 



logA^ 



(59) 



which is the same as the IR divergence in continuum 
QCD, using the same regulator in both theories. In 
the matching coefficients between lattice mNRQCD and 
QCD the logarithmic dependence on the gluon mass will 
cancel out and we can set A = at the end of the calcu- 
lation. 

We discuss three approaches to verify that this same 
divergence is present in the full lattice Feynman integrals. 



1. Infrared subtraction function 

The first approach is to construct a suitable subtrac- 
tion function which can be integrated analytically and 
has the same infrared behavior as the lattice integrand. 



The subtracted lattice integral is then infrared-finite and 
the full result can be obtained by adding the analytical 
expression for the integral over the subtraction function. 
This method was also used in the current matching in 
Ref. [H. 

Only the wavefunction renormalization (in Feynman 
gauge) is infrared divergent. All other renormalization 
constants are IR finite and can be computed directly. 
To construct a suitable subtraction function fi^^^) for 
dZ^i, we start from the continuum integral in heavy quark 
effective theory. (Note that in principle y(''"t>) is arbitrary 
as long as it: agrees with the lattice integrand for small 
loop momenta fc; is ultraviolet-finite in d = 4 dimensions; 
and can be integrated analytically.) 

The logarithmic UV divergence can be regulated with- 
out changing the infrared behavior by replacing 



2jm 



iv ■ k 



(fc- 



(60) 



in the (Euclidean) heavy quark propagator. The result- 
ing integral (which is not restricted to the Brillouin zone) 
is readily evaluated and gives 



6Z. 



(sub) 



Stt 



logA2 + C'(A/m). 



(61) 



This is exactly the logarithmic divergence found in (j59p . 
The subtracted integral SZ^ is evaluated numerically, de- 
fined through 



SZ^ 



(2 



^(eBz(A;)/('^*)(fc)-/(-'^)(fc) 



5Z,i 



Svr 



logA^ 



(62) 



where 0Bz(fc) is equal to 1 inside the Brillouin zone and 
vanishes for any \k^\ > ir/a. 

While this method is easy to carry out for the case 
of the self-energy and vertex correction calculations, it 
becomes increasingly complicated when considering other 
calculations. 



2. Direct calculation for different X 

The alternative, more generic way of isolating the IR 
divergent behavior is to run our integration for different 
values of A and then obtain the desired log A^ behavior by 
numerically fitting a line through the points. For exam- 
ple, in Fig. [6] we show the wavefunction renormalization 
for A^ varying from 10~* to 10~^. Using a logarithmic 
scale on the horizontal axis we see a very clear linear be- 
havior, which demonstrates the desired dependence on 
logA^. The fit to Co + Cm log A^ yields, with a per 
degree-of-freedom of 0.17, Cr = -0.21220(14) log A^, 
which agrees well with the analytic result —2/ (Stt) = 
-0.2122.... 
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FIG. 6: Plot of the wavefunction renormalization for differ- 
ent values of the infrared regulator A. The results exhibit a 
very clear dependence on log A^, and by fitting a straight line 
through the points one obtains Cm = —0.21220(14) ~ — ^ 
and Co = -0.1291(18) with a reduced of 0.17. The value of 
Co agrees well with —0.13124(52), the value of 5Z^ obtained 
directly with A^ = 10~^ (a less precise value is obtained by 
adding fio and Qi in Table |XVII|I . The inset shows the data 
divided by the fit. We use the simple action setup described in 
Appendix IE II The frame velocity is u = 0.3 in this example. 



H. Tadpole improvement 

The tadpole improvement of the action was described 
in Section IIV Bl We define uq to be the mean-link in 

(2) 

Landau gauge. In perturbation theory uq = 1 — UsUq + 

(2) 

. . ., with Uq = 0.750 for the Symanzik-improved gluon 
action [s^. Mean- field corrections are then included as 
counterterms in the action. This leads to 



(tadpole) 



(63) 



where ^7^*'^'^^°'°^ are the resulting tadpole factors which 
we give explicitly below. 

We choose the form of the time derivative in ([^5]) so 
that the wavefunction renormalization is immune from 
mean- field corrections [T^ ]. Thus we expect (and, in- 
deed, find) that the tadpole improvement contributions 
to ilo and fli are exactly equal and opposite. The ap- 
proximate reparametrization invariance implies that the 
radiative corrections to Zp should be small, which sug- 
gests the tadpole corrections to and to H.^ should be 
very similar. Again, we find this to be the case. 

The computation of the tadpole factors was checked in 
two separate calculations. We find 



The latter method can be applied to all kinds of 
calculations such as current matching calculations in 
mNRQCD. It can also be used when the expressions for 
the diagrams are so complex that obtaining the infrared 
counterterms analytically is not feasible. For the inte- 
grands considered here this method is not very resource- 
or time-intensive; even a preliminary investigation, with 
short integration runs and a small number of sampling 
points, can yield a plot with a very good fit, demon- 
strating clear log A^-dependence. For more complicated 
integrands, subtraction functions may still be necessary: 
the computer time required for VEGAS to sufficiently re- 
duce the statistical errors as we lower may well be 
prohibitive and, in addition, strong IR divergences can 
confuse the importance sampling used by VEGAS. 



3. Twisted boundary conditions 



Alternatively infrared divergences can be regulated by 
working on a lattice of finite size and using twisted peri- 
odic boundary conditions 0, which provide a lower 
momentum cutoff. We have successfully implemented 
and tested this method but will not discuss it further 
here. More details will appear in a forthcoming publica- 
tion. 



n, 



(tadpole) 
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(tadpole) 



(tadpole) 



^(tadpole) 
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5 ^^3 
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3-6v^ + 57)'* 
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7^ 777^ 
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1 

477 
1 
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1671^ 



27;^ 
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7777 
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,(2) 



27;2 



V 

6772 



For ?; = these expressions reduce to the ones obtained 
in [57] . Numerical values are given below in Table IIIII 

We give the corresponding expression for an alterna- 
tive treatment of tadpole cancellation in Appendix ID] and 
list the tadpole improvement factors for other, simpler 
actions in Appendix [E] 
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I. Perturbative results 

In this section we present one-loop perturbative re- 
sults for the renormahzation of the mNRQCD propaga- 
tor. Further results for a variety of simpler mNRQCD 
actions are given in Appendix [El 

To obtain agreement with our numerical simulations, 
it is important that we use the Liischer-Weisz gauge ac- 
tion [58l . I59I which is used for the generation of MILC 
lattices 60] . For the heavy quark self energy at one- loop 
level, this action is equivalent to the tree-level Symanzik- 
improved gauge action 

+ 0{as) , 
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^2 




.00 


-2.36685(40) 
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-2.36672(39) 
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-2.26205(38) 
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2.8646(14) 
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2.7636(13) 
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-2.06318(35) 


1.75964(61) 


2.5330(15) 


2.6023(12) 


.60 


-1.91598(33) 


1.62928(62) 


2.3020(17) 


2.4059(12) 


.70 


-1.72666(31) 


1.46150(63) 


2.0220(20) 


2.1623(11) 


.75 


-1.61272(30) 


1.36128(65) 


1.8614(24) 


2.0247(11) 


.80 


-1.48224(28) 


1.24847(69) 


1.6828(29) 


1.8794(11) 


.85 


-1.33083(27) 


1.12528(82) 


1.4925(41) 


1.7275(12) 


.90 


-1.15125(25) 
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1.2930(68) 


1.5972(15) 


.95 


-0.92738(24) 


1.0698(21) 


1.236(19) 


1.6559(25) 



where P, R are 1x1 and 2x1 Wilson loops respectively. 
0{as) denotes possible radiative corrections and tadpole 
improvements of the action that only contribute at higher 
loop orders in the perturbative calculation of the heavy 
quark self energy. 

For the squared gluon mass we choose a value of 
— 10^^. The infrared- finite part of the wavefunc- 
tion renormahzation was extracted using a suitable sub- 
traction function and we also checked that our results 
are indeed infrared- finite by varying A. The stability pa- 
rameter is n = 2 and for the heavy quark mass we use 
TO = 2.8. 

In Table[n]we list numerical results for ilj for a range of 
frame velocities before including mean-field corrections. 
We only give the finite parts of the ilj, the infrared di- 
vergence —2/ (Stt) log A^ is not included in the results for 
r^i, and f2u. 

We give results for the tadpole improvement coef- 
ficients i7(*'^<^P°i«) in Table [m] (see Table IXVll in Ap- 
pendix [D] for an alternative prescription) . Finally we 
show the infrared-finite renormahzation parameters, in- 
cluding mean-field corrections, in Table [TVl and Fig.[7l In 
particular, note that the one-loop coefficient renormaliz- 
ing the momentum is indeed small, as expected from the 
arguments presented in Sec. IV CI 



VI. NUMERICAL SIMULATION RESULTS 

In addition to the perturbative calculations described 
in the previous sections, we have performed a wide 
range of nonperturbative computations with the full 
mNRQCD action on unquenched gluon configurations. 
We have computed two-point correlation functions for 
various heavy-heavy and heavy-light mesons at different 
momenta and boost velocities. These allow the extrac- 
tion of both energies and amplitudes. From the combi- 
nation of simulation energies at different momenta, we 
have obtained nonperturbative results for the external 



TABLE II: Infrared-finite part of fi^ for the full 0{l/m'^ ,vf^i) 
action. The gluon action is Symanzik-improved with = 
10~^ and we use m — 2.8, n = 2. Mean-field corrections are 
not included and the errors shown are purely statistical from 
the VEGAS integration. 
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TABLE III: Tadpole improvement corrections j7(*'''^p°'°) fgj- 
the full 0{l/m^ ,vfgi) action. The heavy quark mass is m = 
2.8 and the stability parameter n — 2. Note that q(''''*p°''=) — 

r-v(tadpolc) 



momentum renormahzation, the energy shift and the ki- 
netic masses of the mesons. We have also examined the 
dependence of several energy splittings on the boost ve- 
locity. In addition to these spectral properties, we stud- 
ied the behavior of decay constants. 

The next section describes the simulations with heavy- 
heavy mesons and is followed by a section on heavy-light 
mesons. All results are given in lattice units. 
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TABLE IV: Heavy quark renormalization parameters for the full 0{l/m^ ,vf^i) action. The gluon action is Symanzik improved 
with = 10"*^ and we use m = 2.8, n = 2. All mean-field corrections are included and the results are infrared-finite. The 
errors shown are purely statistical from the VEGAS integration. 




Name n L S J P C r(r) 
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0.4 0.6 
frame velocity v 

FIG. 7: Heavy quark renormalization parameters for the full 
0{l/m^ ,vfgi) action. The gluon action is Symanzik-improved 
with = 10~^ and we use m — 2.8, n = 2. SZ^ is the 
infrared-finite part of the wavefunction renormalization as de- 
fined in (|62|l . Violation of reparametrization invariance is very 
small, as indicated by the smallness of 5Zp. 



A. Heavy-heavy mesons 



Methods 



We begin by constructing "smeared" interpolating 
fields for quarkonium. To demonstrate the effect of the 
moving NRQCD field redefinition, we start the construc- 
tion with the QCD fields 5*, A meson with momentum 



Vb{lS) 1 - + exp[-\r\/rs]f 

7^b{2S) 2 - + [1- |r|/(2r,)] exp[-lrl/(2r,)] 7^ 

T(15') 10 11 Gxp[-|r|/r,] ^ 

T(25') 2 11 [1- |r|/(2r,)] exp[-|r|/(2r,)]7^' 

X6i(lP) 1 1 1 1 + + cxp[-|r|/(2r,)] (r X 7)Vr, 



TABLE V: Some (continuum) quantum rmmbers and smear- 
ing functions for the bottomonium system. 



p can be obtained from 
Or(p,-r)= ^ l'(a;i,T)r(a;i - a;2)*(a;2,T) 



-tp- 



where r(r) is a Dirac-matrix- valued smearing function. 
We do not include gauge links in r(r); instead we fix 
the gauge configurations to Coulomb gauge. The (con- 
tinuum) quantum numbers and corresponding functions 
r(r) used in the simulations are listed in Table fVl 

We now express ^ and 5* through the tree-level moving 
NRQCD field redefinition. To lowest order one has 

^-(2;) = — 5(A)e-''^'"(-*^-''-=")'^V^(a;), 
Let us, for example, consider the T states with polariza- 
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tion j = 1,2,3. We allow different smearing at source and we obtain 

sink, so that T^^ir) = 7^ fsc{r) and rsk(r) = fsk{r)- 

Using 

W)"/' SiA) = A\r, (64) 



X A\A^ ^^l{xi,T)a^i;^{x2,T)^lix'2,T')a"'^^{x[,T')+ ... 



(65) 



(no summation over j here) where 



k = p — 2'jmv. 



(66) 



The ellipsis in (|65p denotes terms that do not contribute to the connected meson correlator for r > r'. The correlator 
is then given by 



U X\ ,X2:Xi ,xL 



fskixi - X2)e''^ 2 fsc{x[-X2) 



X „Tr I fx' 



G^^;''((a;2,r),(a3^,r')) a™ G^-" ((a^i, r), (a.^, r')) ), (67) 



where we average over TV gauge configurations U. The 
trace is over color and spin indices. We have also used 
equation ((25|) to express the antiquark green function 
G^'" in terms of the quark green function G^'~'" with 
the opposite boost velocity. 

The summations over all quark and antiquark source 



locations would render the lattice computation too ex- 
pensive. Therefore, using translation invariance, we re- 
move the summation over the antiquark source location 
x[. Furthermore, we remove the factor of Q-'^f'^i'^-'^ ) 
which corresponds to the tree-level energy shift. Hence, 
the quantity 



G(r.k,r.e,fc,T,T') = E e-'"^ fsk{x^ - X2) 



xA^A\,Tr(a' [g^;'' ((a;^, r), (a.;, r')) 



<r ((^i,r),(a3;,r')) 



with 



Gl^: {{X2,r), ix[,r')) = ^ e^^^/.c(a.'i - a=^)G^f {(x2,t), [x'^^)) 



(68) 



is computed on the lattice. The correlator can be 
computed by using the function 

e*'=^/sc(a='i-4) (69) 

as the initial condition in the mNRQCD evolution equa- 
tion (j28p . The momentum-dependent phase factor 
exp(zfc(a;j + x'2)/2) at the source improves the overlap 



with the momentum considered. However, since there is 
no sum over a;^, one may omit this factor to allow the 
calculation of correlators with different momenta from 
the same source. 

In order to maintain the periodic boundary conditions, 
we set /(r) to zero for \r\ > Rg with some cut-off radius 
Rs smaller than half the length of the lattice. 
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On the finite volume lattice with periodic boundary 
conditions, the momentum k takes on discrete values, 
kj — 2tt rij I Lj where Lj are the spatial extents of the 
lattice. However, the physical meson momentum p is 
expected to deviate from the tree-level relation 166]) , since 
mass and velocity are renormalized. One has 



p = 2ZpPo + k with Pq = ^mv. 



(70) 



We fit a matrix of correlators with different smearings at 
source and sink with the functional form 



C(r,k,r,e,fc,r,r')->A^'^(A^ 



-E{r- 



+ 5sk(5SC)*g-(B+ABi + ...+Ai;„)(r-r' 



(71) 



where E is the energy of the meson ground state, A^^ 
and A^"^ are the (real) ground state amplitudes of the 
operators at source and sink and B^, are (real) am- 
plitudes for the n-th excited state, relative to the ground 
state amplitude. We use the constrained fitting method 
described in |6l| , and increase the number of exponentials 
until the fit results and error estimates become indepen- 
dent of ricxp- 

The full (physical) energy differs from the energy E = 
Ey{k) obtained from the fit by twice the mNRQCD en- 
ergy shift, 

Epi,ys = E,{k) + 2a. (72) 

In perturbation theory, one has 

Cy = ZrnZ^-fm + Eq. (73) 

Given expression (|70p for the full (physical) momentum, 
we expect that, up to lattice artifacts. 



E, 



phys 



= J{2ZpPo + k) 



M2 



(74) 



where Afkin is the kinetic mass of the meson. 

Using (|74p . we can obtain nonperturbative results for 
Cy, Zp and Mkin from the energies at various non-zero 
lattice momenta in combination with the energy at A; = 0: 



Cy 



lkl~{E^^{k^)-E^M) 



2 2{Ey{k^) - Ey{0)) ' 
i?2(fc||)-i?2(_fc[|)+4a(i?.(fc||)-i?„(~fc||)) 



(75) 



4fc|| • 2Po 



Mkin = \ iEy{k) + 2CyY - {2ZpPo + fc)2. 



(76) 
(77) 



Here, fc|| is parallel to v, and k± is perpendicular to v. In 
order to fully take into account correlations in the ener- 
gies at different momenta, we use the bootstrap method, 



performing fits on 500 bootstrap ensembles and comput- 
ing the final quantity 500 times. The errors are then 
estimated as the 68% width of the resulting distribution. 

Ultimately we will be interested in semileptonic B de- 
cay matrix elements. As a simpler test we first study the 
decay of the r]i,{lS) meson via a fictitious axial vector 
current. The corresponding decay constant is defined by 



|J^(0)|77f,(15),p) 



(78) 



Here, Jg is the mNRQCD field operator associated with 
the axial current 



'^{x) ^^"Y ^{x). 



(79) 



For simplicity, we have only considered the temporal 
component and, as above, used only the leading-order 
tree-level mNRQCD field redefinition to construct the 
lattice current. To extract the amplitude, we compute 
2x2 matrix correlators with the local smearing function 



r(r) = 5{r) 



(80) 



for the temporal axial current, and the 'qi,{\S) smearing 
function from Table |Vl The product of the ground state 
amplitudes in (j7ip is given by 



1 



-(r;,(15),p|Or.,(0)|0) 



X (0|Or..(0)|ryb(15),p), (81) 

as can be seen from the spectral decomposition of the 
two-point correlator. Using ((78)) with = i?phys, ([TS]) 
and (|8ip. we obtain 



Ey{k)+2Cy 



(82) 



where A = yj^sk/sc |-j^g amplitude from the fit corre- 
sponding to Tsk/sc = '5(r) 7^7*'. 

2. Lattice parameters 

The computations were performed using 400 MILC 
gauge configurations (fixed to Coulomb gauge) of size 
20'^ X 64 with 2-1-1 flavors of rooted staggered light 
quarks, at (3 — 6.76 [60]. The light quark masses were 
my — rud — 0.007 and = 0.05 (in the MILC conven- 
tion for lattice masses). The Landau gauge mean link, 
used in the mNRQCD action, was uq = 0.836. The in- 
verse lattice spacing of these "coarse" MILC configura- 
tions is known to be approximately 1.6 GeV (6^ . 

Heavy quark propagators were computed using full 
mNRQCD lattice action described in section HVl and used 
in the perturbative calculation. The bare heavy quark 
mass was set to m = 2.8, which gave the correct T kinetic 
masses using non- moving NRQCD [g^]- The boost veloc- 
ity was always pointing in the a;-direction, v = (7J,0,0). 
The stability parameter was set to n = 2. 



22 



In order to increase statistics, between 16 and 120 cor- 
relators with different origins {x[,t') spread over the lat- 
tice were calculated and averaged over on each gauge con- 
figuration. These origins were also shifted randomly to 
reduce autocorrelations. The smearing parameter was 
set to 1 for the S wave states and 0.5 for the P wave 
states. 



3. Results 

Results for the r]b{lS) kinetic mass Mkin and the renor- 
malization parameters Zp, Cv are shown in Table IVII 
The energies were obtained from 6-exponential fits to 
2x2 matrix correlators with the r]i,{lS) smearing and 
the local axial current. For the calculation of Cy using 
([75)) . we averaged the results over the 4 different perpen- 
dicular lattice momenta 

fc^e | — (G,±1,0),— (0,0,±l)j. (83) 

The momentum parallel to the boost velocity in (|76p was 
chosen to be fc|| = ^(1,0,0), and in (|77p . for the mea- 
surement of Mkin, we use fc = 0. 

Because the lattice is of finite extent, L = 20 in our 
test case, the estimates for Cy and Zp will be affected by 
the choice of momenta in ([75)1 and ((76|) since the formu- 
lae are accurate only in the limit that the momenta are 
infinitesimal. Note that the uncertainty due to using non- 
infinitesimal momenta will decrease for larger lattices for 
which smaller momenta are available. 

To estimate the size of the resulting systematic error we 
also performed the calculations with the larger momenta 

fc^e | — (0,±2,0),—(0,0,±2)|,fc|| = —(2,0,0). 

(84) 

For Cy, the results from \k±] — 2-k/L agree with those 
obtained from |A;^| ~ At: / L within statistical errors, indi- 
cating that the systematic error is small and does not in- 
crease significantly when increasing the momentum per- 
pendicular to V in the measurement. For the measure- 
ment of Zp at \v\ — 0.6 we find a 6% {2a) change in 
Zp when going from |fc||| — 2'k/L to |fc||| = Att/L. At 
\v\ = 0.4 and smaller boost velocities the results are equal 
within statistical errors. For the kinetic mass, which de- 
pends on both Cy and Zp, we again find agreement within 
statistical errors between the results from the two differ- 
ent momenta for all boost velocities considered. At small 
velocities, we find that both Zp and C^/i'^m) are close 
to their tree-level value of 1, demonstrating that renor- 
malizations are indeed small. 

We also obtained the amplitude for the axial current 
and extracted the pseudoscalar decay constant from the 
same 2x2 matrix fits using ([82]) . For the energy shift 
Cy in ([5^ we used the result from |fc^| = 2tt/L. The 
meson momentum is given by p = 2Zpjmv + k. In the 
following we compare two methods of reaching large \p\. 



■ NRQCD 
« mNRQCD 

0.8 - 

I 

0.6 - 



0.4 - 



0.2 - 



FIG. 8: Heavy-heavy decay constant in NRQCD and 
mNRQCD for different values of the meson's momentum, 
\p\/{2tv/L) = 0...10 (NRQCD) and p ^ Zp 2'ymv for 
|wi = 0.2,0.4,0.6 (mNRQCD). The horizontal line indicates 
the value at p = 0. 



First, at D = 0, i.e. with standard NRQCD, we computed 
the decay constant at large non-zero lattice momentum fc; 
the results are shown in Table IVTIl Second, we computed 
the decay constant with fe = and three different boost 
velocities v; the results are shown in Table IVini In this 
case the uncertainty in Zp leads to an uncertainty in the 
meson momentum. 

A plot of the decay constant against the total momen- 
tum (with Zp from ([7S|) with |fc||| = 2tt/L) for the two 
methods is shown in Fig. [S] The decay constant is a 
Lorentz scalar and should be independent of the momen- 
tum. However, with NRQCD we see large deviations due 
to both relativistic and discretization errors. With mov- 
ing NRQCD the deviation is very small, giving evidence 
that the formalism works very well. Small deviations are 
still expected here, since only the leading-order current 
was used; i.e. Tfwt and A^^ were set to unity in (j2ip for 
this calculation. 

Next, we studied the velocity-dependence of various 
energy splittings between the bottomonium states listed 
in Table |Vl For the T and r]b states, we used 6- 
exponential 2x2 matrix fits with the IS" and 25* smear- 
ings; for the Xbi states a 6-exponential single-correlator 
fit with the IP smearing at both source and sink was 
used. The results for the T (25') - T ( 15) , Xfci ( 1 ^) - T ( 1 S") 
and T (15) - 776(15) splittings are hsted in Tables Ull H 
and IXIl respectively. 

Note that the energy splittings are not Lorentz scalars. 
Using ((74|) . we expect that the splitting between two 
states A and B at zero lattice momentum is given by 

E^{0)^E^{0) = ^{2Zp^mvy + {M^J^ 

-^{2Zp^mvY + {Mi^Y. 
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\k±\ = fc|| = 2n/L 


fc 1 1 = Ifeiil = 47r/L 


\v\ 


Zp Mkin Cv/i^m) 


Zp Afkin C„/(7m) 





— 5.974(48) 1.0182(86) 


— 5.979(37) 1.0190(65) 


0.2 


1.008(19) 5.95(10) 1.015(18) 


1.009(12) 5.969(62) 1.017(11) 


0.4 


0.9953(78) 5.931(44) 1.0084(77) 


0.9830(65) 5.954(40) 1.0101(70) 


0.6 


0.898(27) 6.22(18) 1.010(28) 


0.843(27) 6.37(15) 1.011(21) 



TABLE VI: Nonperturbative results (using the 77f,(lS')) for Mkin, Z^^ C„. 



|p|L/(2^) 




/ 








0.4724(23) 


1 


0.31416 


0.4731(23) 


2 


0.62832 


0.4755(24) 


3 


0.94248 


0.4772(43) 


4 


1.25664 


0.4835(77) 


5 


1.57080 


0.4971(78) 


6 


1.88496 


0.5209(46) 


7 


2.19911 


0.5527(44) 


8 


2.51327 


0.6006(45) 


9 


2.82743 


0.6740(49) 


10 


3.14159 


0.715(29) 



TABLE VIL r?6(lS') decay constant with standard NRQCD 
(i.e. = 0) computed with several values of meson momen- 
tum IpI by varying jfe|. 



\v\ 



A£.„(0) 



AJS^(O) 
A£;o(0) 



0.0 0.2703(89) 1 

0.2 0.264(12) 0.976(56) 

0.4 0.270(23) 0.998(91) 

0.6 0.227(57) 0.84(21) 

TABLE X: Xb\iXP) - T(IS') energy splitting as a function of 
the boost velocity. 



A£„(0) 



Ag„(0) 
ASo(O) 



0.0 0.031469(98) 

0.2 0.03039(20) 

0.4 0.02837(85) 
0.6 0.0281(28) 



0.9656(71) 
0.901(27) 
0.894(88) 



/ 



0.4724(23) 

0.2 1.152(22) 0.4739(38) 

0.4 2.433(19) 0.4810(36) 

0.6 3.77(11) 0.499(11) 

TABLE VIIL 776(15') decay constant with mNRQCD at fc = 
computed with several values of meson momentum \p\ by 
varying \v\. 



\v\ 



A£,(0) 



Ag^O) 
Ai;o(0) 



0.0 0.3334(68) 1 

0.2 0.329(10) 0.986(37) 

0.4 0.320(15) 0.958(48) 

0.6 0.20(11) 0.59(33) 



TABLE IX: T(2S) 
the boost velocity. 



T(IS') energy splitting as a function of 



TABLE XI: T(IS') 
the boost velocity. 



ryb(lS') energy splitting as a function of 



Uj 



Uj 





' 1 ' 


1 










■ Y(2S)-Y(1S) 


1 






, X„{1P)-Y(1S) 








« Y(1S)-riJ1S) 








— 1 -0.5v^ 








, 1 , 


1 



FIG. 9: Bottomonium energy splittings relative to « = as a 
function of the boost velocity. Points are offset horizontally 
for legibility. The data agree with an estimate for the leading 
dependence (see text). 
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\v\ 


A TT' /'r\\ 1 

A_£;„(0)|i 


a_e;„(0)|2 


A 7T1 /'r\\ 1 

AEi,(0)\3 





-0.000009(63) 


-0.000039(68) 


0.000053(73) 


0.2 


-0.00012(26) 


-0.00005(28) 


0.00017(30) 


0.4 


-0.00046(56) 


0.00055(62) 


-0.00010(57) 


0.6 


-0.0176(96) 


0.0107(62) 


0.0069(75) 



TABLE XII: Dependence of the T(15') energy on the polar- 
ization direction. Ai5„(0)|-, is the difference between £'„(0)jj 
and the polarization-averaged energy. 



If we set Zp = I and expand the splitting at velocity 



V relative to v 
obtain 



in powers of the boost velocity, we 



2m^ 



A,/ -4 71/rB 



that is, we expect a quadratic decrease like 1 — 0.5|t>p. 
The numerical results, shown in Fig. [51 are consistent 
with this estimate as desired. 

Finally, for the T{1S) meson, we studied the depen- 
dence of the energy on the polarization direction. If 
moving NRQCD works well, then there should be no 
difference for polarizations parallel and perpendicular to 
the boost velocity. In Table IXIII we show the difference 
between the energy with definite polarization direction, 
Ey{0)\j and the polarization-direction-averaged energy 
i(£;^(0)|i + -B„(0)|2 + E^{0)\3). No significant depen- 
dence on the polarization direction can be seen (except 
maybe at w = 0.6, where a l.Scr deviation in the energies 
was found). 



B. Heavy-light mesons 

1. Methods 

Starting with the standard Dirac fields, we construct 
interpolating fields for the Bg and B* mesons with mo- 
mentum p from 



OriP,T) 



^i{x, T)r{x - y)*H(y, T)e-'P y , (85) 



where '3// is the Dirac spinor for the valence strange 
quark and is the Dirac spinor for the b quark. We 
use r(r) = 7^ /(r) for the Bg pseudoscalar meson, 
r(r) = 7^ /(r) with j = 1,2,3 for the B* vector meson 
and r(r) = 7^7° /(r) for the computation of the decay 
constant /s^ . We compute 2x2 matrix correlators with 
Gaussian and local smearing, f{r) — e^l'"' , S{r). 
In terms of the standard Dirac propagators, the two- 



point function reads 

{OrAP.r)Oljp.T'))^^Y. E 



x,y,x ,y 



X e 



Tr [ T,^{x - y) Gi {x' , x) TI{x' - y') Gh {y, y') ] 

(86) 



with X = (a;,T), y = (y,r), x' = {x',t'), y' = {y',T'). 
For T > t' , the tree- level leading-order mNRQCD field 
redefinition pT|) leads to the following expression for the 
b propagator: 



Gh (y, y') = -e 



1 

— ( 

7 



-7m(r— r') +17771 V - (y — y') 



,^(^)(^G,„(,, y') Oj^^^^_ 

For the light quark, we use the ASQTAD staggered 
fermion action 34]. The 4-component naive light quark 
propagator can be obtained from the 1-component stag- 
gered propagator G^{x' ,x) via 



with 



Gi{x\x)= G^{x\x)®n{x')9.\x) 



(87) 



(Recall our convention for the Dirac matrices is as given 
in Appendix [XI) We also employ 7^-hermiticity 



Gi{x',x)^fGj{x,x')f 



(89) 



to interchange the points x and x' for the light quark 
propagator. As before, we remove the factor of 
g-7m(r-r ) |.]^g Summation over x'. 

In the case where Tsk and Fgc contain the same Dirac 
matrix, we arrive at the following expression: 

C(r,k,r,c,fc,r,r') = lj2^J2f,^(x-y)e-^'^y^ix,x') 



X Tr 



G\{x,x') s{m{^')^H^)s{i^)r^^^l ""'^ I 



(90) 

with k = p ^mv and 

(y, = E /(^' - y')e^'^-«'G^„ (y, y'). 

y' 

The phase factor ri{x,x') in ([90)1 depends on the Dirac 
matrix in Fgi^ and Fgc. It is given by 

1 for 7^, 

r]{x,x') = { (-l)^3-^j for 7^', 

(_l)Ej(^i+4) for j5^o^ 
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As before, we set /(r) to zero for \r\ > Rg with some cut- 
off radius Rg smaller than half the length of the lattice. 

The staggered/naive light quark action used here suf- 
fers from the doubling problem. As shown in [s^l, the 
spatial doublers do not contribute to the correlators. 
However, the temporal doubler leads to a coupling to 
additional opposite parity states, which manifest them- 
selves as oscillating exponentials in the correlators. We 
therefore fit the heavy-light correlators to 



-E(t-t') 



.+AE„){t-t') 



+ {-iy-^'+^A''^{A"') 



^sk^^scy^-{E+AEi + ...+AEr„)iT-T') 



m—1 



The quantities C„, Zp ,Mkin and the decay constants /s, 
/b^ can be extracted in a completely analogous manner 
as for the heavy- heavy- mesons, with the replacements 
2Cv Cv and 2Po — > Pq, since now there is only one 
heavy quark. 



2. Lattice parameters 

The heavy-light simulations have been performed with 
the same gauge configurations as the heavy-heavy simu- 
lations, and the same heavy-quark action and parameters 
were used. Again, the boost velocity was always pointing 
in x-direction, v = (t),0, 0). The valence strange quark 
mass for the Bs and B* mesons was set to 0.040. Four 
staggered propagators with source times r' — 0, 16, 32, 48 
were used for each gauge configuration. Both forward- 
and backward-propagating meson correlators were com- 
puted to increase statistics. The smearing parameter 
was set to 2.5. 



V 


7 




I \ }'''') 







3.37(15) 


1.002(52) 


0.2 


1.05(15) 


3.72(47) 


1.13(16) 


0.4 


1.05(18) 


3.66(68) 


1.10(23) 



TABLE XIIL Bs results for Mkin, Zp, Cv. 
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local - smeared 
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FIG. 10: Bs matrix correlators at fe = and v = (upper 
panel), |?;| — 0.4 (lower panel). 



3. Results 

Results for the B^ kinetic mass Mkin and the renormal- 
ization parameters Zp, are shown in Table IXIIII The 
energies and the amplitude required for the calculation 
of the decay constant were obtained from 8-exponential 
(4 of which are oscillating) fits to 2 x 2 matrix correla- 
tors with the Gaussian smearing and the local axial cur- 
rent. Two sample plots of these correlators at d = and 
V = 0.4 are shown in Fig. [TOl This also demonstrates 
the worsening of the signal-to-noise ratio as the boost 
velocity increases, in accordance with (U). 

For the calculation of C„ , we again averaged the results 



over the 4 different lattice momenta perpendicular to v 
fci = ^(0,±l,0), ^(0,0,±1), (91) 

and the momentum parallel to the boost velocity re- 
quired for the determination of Zp was chosen to be 
fell =^(1,0,0). 

As expected, the statistical errors are larger than for 
the heavy-heavy mesons, partly due to a much smaller 
number of origins (four) per gauge configuration. The 
results for Zp and agree with those obtained using 
heavy- heavy mesons in section IVI A 31 

The results for the decay constant fs^ at fe = and v = 
0, 0.2, 0.4, 0.6 are hsted in Table IXTVl and plotted against 
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/(fc = 0) 



0.1626(27) 

0.2 0.576(11) 0.1608(52) 

0.4 1.2163(96) 0.1634(94) 

0.6 1.885(57) 0.174(17) 

TABLE XIV: Bs decay constant (unrenormalized, and in lat- 
tice units) with mNRQCD at fe = 0. 
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FIG. 12: Renormalization Zp of the external momentum 
Po ~ ymv. We show perturbative and nonperturbative re- 
sults for heavy-heavy and heavy- light mesons, with a slight 
horizontal offset for legibility. The uncertainties shown on the 
data points for the perturbative results are purely statistical 
due to the VEGAS integration. The strong coupling constant 
is taken to be Os = av{q*) and the error band is obtained by 
varying the matching point in the range [q* /2, 2q*]. 



VII. COMPARISON OF PERTURBATIVE AND 
NONPERTURBATIVE RESULTS 



FIG. 11: The Bs decay constant at fc = and v — 
0, 0.2, 0.4, 0.6 plotted against the total momentum p — 
Zpymv + k. The horizontal line indicates the value at u = 0. 



the total momentum in Fig. [TTJ In the calculation of the 
decay constant, we used and Zp determined from the 
776(15) dispersion relation since this is more precise. We 
find that the decay constant is independent of the boost 
velocity within statistical errors. (Even when working 
with non- moving NRQCD, the discretization errors in 
the heavy-light decay constant do not appear to grow 
as severely with momentum (63j as in the heavy-heavy 
decay constant (Fig.[H]).) 

We also computed the B* — Bg energy splitting as a 
function of u; the results are shown in Table IXVI The 
statistical errors are so large that no definite statement 
can be made about the velocity dependence. 



V 


A£;„(0) 


A£40) 


0.0 


0.0261(35) 


1 


0.2 


0.0262(65) 


1.00(28) 


0.4 


0.0310(80) 


1.18(34) 



TABLE XV: B* — Bs energy splitting as a function of v. 



In the following we compare our perturbative results 
given in Section IV II to the nonperturbative numbers ob- 
tained in Sections IVTXgl and rvTB3] 

We use the str ong coupling constant defined in the 
potential scheme [36| and choose g* (for each quantity 
and each value of v) using the Brodsky-Lepage-Mackenzie 
procedure [63 |. The q* values range approximately be- 
tween 0.5/a and 3/a. As a reference. 2/a = 3.2 GeV 
on the coarse MILC configurations [64|. Using the run- 
ning of the strong coupling constant ay (q) [65j this gives 
ay(2/a)«0.3. 

In Figs . [T^ and [T5I we show both perturbative and non- 
perturbative results for the renormalization of the exter- 
nal momentum and the energy shift between QCD and 
mNRQCD (see Section IV Bl) . The discrepancies we find 
at V — 0.6 indicate sizable higher order loop contribu- 
tions as V grows. High-/? simulations verify the one-loop 
perturbative calculation as described earlier, and prelim- 
inary estimates of the gluonic (i.e., quenched) two- loop 
contribution using high-/? simulations show that higher- 
order loop corrections reduce this discrepancy; further 
work is in progress and will be presented in a forthcom- 
ing publication. 



VIII. CONCLUSION 

We have derived the mNRQCD action through 
O (1 / , v^^i) and discretized it with errors starting at 
0{asa^) (tree- level errors begin at 0{a^)). The one-loop 
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FIG. 13: Renormalization of the energy shift Cv compared 
to the tree level value ym. We show perturbative and non- 
perturbative results for heavy-heavy and heavy-light mesons, 
with a slight horizontal offset for legibility. Uncertainties are 
presented as in Fig. [I2l 



renormalizations of the wavefunction, the external mo- 
mentum, the frame velocity, and the energy shift Eq have 
been computed and presented here. In the cases of the 
external momentum and the energy shift, we compared 
perturbative and nonperturbative results. Nonperturba- 
tive calculations of heavy-heavy meson and heavy-light 
meson properties were undertaken, with the aim of test- 
ing the specific action and the general method. Fig. [5] 
is particularly instructive; it shows the reduction in dis- 
cretization errors obtained by using mNRQCD compared 
to non-moving NRQCD to compute the fictitious jyf, de- 
cay constant. Whether mNRQCD will prove indispens- 
able in determinations of heavy-to-light form factors is 
still to be seen. Nevertheless, lattice calculations of these 
form factors are a pressing need, and the more tools we 
have at our disposal, the more quickly can we understand 
and reduce the errors in our calculations. In particular, 
these methods will enable us to explore the ^ limit 
needed for the rare decay B —^ K*j while maintaining 
control over lattice discretization errors for the light vec- 
tor meson. In future work we will employ mNRQCD, 
and other tools, to move towards this goal. 
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APPENDIX A: NOTATION 

In this Appendix we summarize for convenience our 
choices of notation and convention detailed throughout 
the main text. 



• Lorentz boost: 



A 



7 



7 v-^ 

with 7 = (1 - t;2)-i/2 
• gamma matrices: 



7 V 

1+7 



7 



7 



-a" 



*7°7^7^7"^ 




with the Pauli matrices cr-'. We define a'^ = l2x2- 
• spinorial Lorentz boost: 



^(A) 



1 



/ 1+7 J cr ■ V 
v/2(l+7) \1(T V 1+7 



• covariant derivatives and field strength tensor: 

d_ 

9a; 



^1^ = ^+ '9A^ 



• chromoelectric and chromomagnetic fields in 
Minkowski space: 



Ek — Fok, 



• chromoelectric and chromomagnetic fields in Eu- 
clidean space: 



Ek 



4k, 



-7:^jkiFki 
2 
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APPENDIX B: REMOVING TIME DERIVATIVES 
IN H AT ORDER 1/m^ 

In this section we show in detail how additional time 
derivatives can be removed from the mNRQCD La- 
grangian at 0(1/™^). In particular we give an explicit 
expression for the operator V in (fT7|) . 

The field redefinition ([T7]) results in 



+ 0{l/m' 



On 



1 

7m 



o 



1 



(2)1 



(7m)' 



7O 



(2)2 



(2) 



with 



0(2)1 = 



0(2)2 = 0(i)2+{F, Oo}, 

and we need to write 0(i)2 = 0(2)2 — {V, Oq} with some 
operator V such that 0(2)2 does not contain time deriva- 
tives. We will treat the different terms in 0(i)2 (see ([T5)) ) 
individually. Note that the last term, — { — ^C/^, Oo}, is 
already in the desired form. The time-derivative in the 
original O2, defined after (fH|) . can be treated as follows: 

|7 7%MS^A",{i^o, 



(Bl) 



Next, using 



we obtain 



U = 



\h'-l)Oo+'-j°v-D 



(B2) 



UOoU - i{C/2, Oo} + ^[U, [Oo, U]] 



1 



{U', Oo} 



{U^ Oo} 



Oo, i(7'-l)Oo + ^7%-£> 



if{Do + v-D), - fv-D 



1 



= ^{U'.Oo] 
- \ {U\Oo] 



[U, [Do,v-D]] 



16 



1, 



-U\ Oo 



(B3) 



and 

{U, Oi} = 



J(7'-l)Oo + ^7%-r>, Oi 



^7%.D, Oil- 1-1(7^-1)01, Oo 



^7%.D, 0(i)i[- U7°t^-A {[/, Oo} 



1 



-(7^-1)01, Oo 



-j^'v-D, 0(i)i 



-r ''^•£>, Oo 



^7%.r>, [/Ui(y-i)Oi, Oo 



(B4) 



Let us now consider the nested commutator in (IB4I 



i7, 



-7%-D, Oo 



-7° i7°Do 



^7°((7^-l)i?o + (7^ + l)^-£'),ft^-£; 

= -§7° ((7' - i)i?o' + (7' + 1)^-£>^') (^ • £;) . 

We conclude from jHI]), (|B3]) and (jB4| that 
y - -f7e,feiSJA\i?,'+(^7%.r>, [/ 



1 



(72 - l)Oi - C/2 



(B5) 



and 



O, 



(2)2 



-ze,ki^^A\{vD, El} 



16 ' 

i 



-lljfT, -{DxE' ^E' xD) 



^°{v-D, S- (v X £;')} 



8(1 + 7) 



16 



7" (Dg^ -t;-£>^") (t;-£;). 
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APPENDIX C: LATTICE DERIVATIVES AND 
FIELD STRENGTH 



In this section we give explicit expressions for the dis- 
cretized derivatives we use in our lattice action, Eqs. pop . 
pip . All expressions are constructed from the elementary 
forward, backward and symmetric derivatives 



A+ij{x) = Uf,{x)il;{x + fi)-'ijj{x), 
^'ili'i^) = ^{x) -U-f,{x)il){x - (i), 

= \ [U^{x)i^{x + fL)-U-^{x)^{x-ii)]. 



For performance reasons, we construct higher-order oper- 
ators to be maximally local by balancing the occurrence 
of these three types. We also symmetrize the expressions. 

Unimproved derivatives: 



A(2) 
A(2) 
A(3) 
A(4) 



\ «^"^;S'(A+A±Ar + ATA±A+) 



3 



\ ^ t;^Vt;'z;'"(A+A-A+A; 



j.k,l,m—l 



+ATA+ArA+) 



Improved derivatives: 



A± = A± - -A+A±A- 

J 3 Q 3 J 3 



A(2) ^ A(2) - A Va+AtA+A~ 

12 ^ 3 3 2 3 

A(^) = A(^) + i^.^«'=A;A7A+A- 

V vH'' fA+A^A+Ar + ATA+ATA+ 

Yl ^ 3 3 3 k ^ ^] ^3 3 k 

3,k=l 

+A;A,-A+A, +ATA+A,A+) 
Unimproved adjoint derivative: 

U^{x)Fp„{x + fL)Ul{x) 



U.^{x)Fp„{x ^ fi)UlJx) 



Improved field strength tensor: 

1 



U^{x)F^Ax+fi)Ul{x) 
+ U-p{x)Fp,{x-fi)Uip{x) 
- (m ^ 



where 
Ff^Ax) 

with 



2^ 



'n^Ax)-nl,{x)), 



\ Ua{x)Ui3{x+a)U.a{x+a+f3)U.i3{x+f3) 



APPENDIX D: TADPOLE IMPROVEMENT 



In the perturbative calculation it is possible to explic- 
itly work out every path appearing in the evolution and 
cancelling the tadpole factors which appear in every in- 
stance of U^{x)Uj^{x). Here we give analytical expres- 
sions of the tadpole improvement corrections for this case 
for the fuU ©(l/m^, w^^J action. 

Numerical results for m = 2.8 and n = 2 can be found 
in Table IXVII and should be compared to Table IIIII 



(tadpole) 




= -VL\ 



(tadpole) 



V 



3 768 1024 



2688- 852i;2 + 11t;4_i3„6 





7687m 






3456- 


- 4920i;2 -)- 24971;-* - 


- 264wf^ - 


f 15w*^ 




307272m2 






516- 


1264z;2 -1- 1058?;4 + 


275^6 - 


15?;* 




7687'^to3 






-591 


+ 1460u2 - 135817" 


+ 4481-*^ 


+ 5w® 



25674m4 

81 - 216t;2 + 246w* - I2?>v^ + 2bv^ 
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V 


A(tadpolc) , (2) 


^(tadpole) /^(2) 




0.00 


2.10610 


—3.14336 




0.01 


2.10600 


—3.14319 


—3.14321 


0.10 


2.09592 


—3.12639 


o 1 o or\o 

—3.12898 


0.20 


2.06507 


—3.07557 


o r\C} f TO 

—3.08573 


0.30 


2.01267 


—2.99112 


o r\ 1 o o 1 

—3.01321 


0.40 


1 no TOO 


—2.87355 


o n 1 no o 

—2.91083 


0.50 


1.8J666 


o TO /I n /I 
— z.7z4U4 


—2.77787 


0.60 


1.70836 


— 2.54438 


—2.61351 


0.70 


1.54900 


O O O T /I O 

—2.33743 


—2.41672 


0.75 


1.45625 


-2.22476 


-2.30557 


0.80 


1.35355 


-2.10633 


-2.18513 


0.85 


1.23910 


-1.98182 


-2.05402 


0.90 


1.10911 


-1.84923 


-1.90880 


0.95 


0.95262 


-1.70037 


-1.73901 



(2) 



TABLE XVI: Tadpole improvement corrections {^(''"^p"''^) foj- 
the full 0(1/771^,11*31) action and cancellation of U^Uj^ as de- 
scribed in the main text. The heavy quark mass is m = 2.8 



and the stability parameter n = 2. Note that Q,\ 

^(tadpole) 



(tadpole) 



-f2, 



that the difference is of the order of 10% in the one-loop 
coefficient, see Tables IIIII and IXVII We conclude that it 
is sufficient to avoid multiplying Ufjjj^ by I/uq within 
and 5H separately. 



APPENDIX E: FURTHER PERTURBATIVE 
RESULTS 



In this appendix, we present one-loop perturbative re- 
sults for the renormalization of the mNRQCD propagator 
for various simpler forms of the mNRQCD action. 



1. Simplest case 

We considered the simplest, unimproved mNRQCD ac- 
tion, i.e. 







Hq = -iv ■ - 
SH = 



± A(2) - A?) 



27TO 



(tadpole) 



(2)/ 5 7«2 
"V 3+32 



512 2048 



-10880 + 4480i;2 - 215?;'' 


+ 35v'^ 




30727m 






-12480 + 10288w2+4321i;4 


-360^6^ 


-f ISw^ 


61447^771^ 






2412 - 4864i;2 + 3974i;4 -f- 


311i;6 - 




15367^771-^ 






-879 + 2100w2 - 1982w4 -f 


640i;^ 4 


- 5v^ 



512747774 

81 - 216i;2 + 246i;4 _ 128w^ + 25v^ ^ 



1287^7^ 



5 llw^ 29^4 v'^ 



^(tadpole) ^ ,(2)/ ^ + + + 

° V 3 48 1536 



2048 

-5440 + 1860i;2 - 51i;4 -|- I6v^ 
153677)7 

12480 + 712v^ - 352l7;4 + 320v^ - 15i;^ 



61447^7(72 






2412 - 3016^2 + 2306^4 4 


- 299v^ - 


157;** 


15367^777^ 






-879 + 18127;2 - 16147;4 - 


f 5447;'' -1 


- 5v^ 



81 - 2167;2 + 246i;4 - 1287;^ + 25v^ - 



1287^7775 J 

It should be noted that the expressions for partial can- 
cellation are significantly simpler. Numerically we find 



coupled to the Wilson gluon action. The gluon propaga- 
tor in Feynman gauge is 



D-i(fc) = 4^sin2^ + A2 



2 - 7« - w^^ + 4 X! ™^ 17 + 



with w — e^^^ . The gluon mass was set to = 10 ^. 

The case 5H = is very simple, as all propagators and 
vertices are diagonal in spinor and color space, and the 
calculations can be performed in reasonable time on a 
workstation. We used a heavy quark mass of m = 2.8 
and the stability parameter is 77 = 2. 

In Table IXVlII we list J7j for this action before includ- 
ing mean-field corrections. We only give the finite parts 
of the fij", the infrared divergence — 2/(37r) log is not 
included in the results for O2 and O^. 

The mean-field corrections, cancelling U^Uj^ factors as 
described in the main text are 



(tadpole) 



_^(tadpolc) ^ ^(2) ( ^ 



3-7;^ 
7777 



(tadpole) 



277 - 1 3 - i;2 



(El) 



277 7777 

whereas the corresponding expressions for tadpole can- 
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V Qo Qi SI2 f2 




0.00 


-2.9851(24) 


2.8619(24) 


3.9967(29) 


— 


0.01 


-2.9879(24) 


2.8645(24) 


3.9987(29) 


4.003(23) 


0.10 


-2.9721(24) 


2.8483(25) 


3.9889(29) 


3.9741(39) 


0.20 


-2.9299(23) 


2.8033(24) 


3.9567(29) 


3.9474(31) 


0.30 


-2.8564(23) 


2.7252(24) 


3.9022(29) 


3.8826(29) 


0.40 


-2.7490(22) 


2.6092(23) 


3.8218(29) 


3.7898(28) 


0.50 


-2.6085(22) 


2.4540(22) 


3.7104(30) 


3.6702(27) 


0.60 


-2.4260(20) 


2.2462(21) 


3.5651(33) 


3.5087(27) 


0.70 


-2.2057(18) 


1.9859(20) 


3.3833(39) 


3.3157(25) 


0.75 


-2.0832(18) 


1.8335(20) 


3.2742(45) 


3.2110(26) 


0.80 


-1.9371(17) 


1.6482(19) 


3.1333(57) 


3.0851(26) 


0.85 


-1.7790(16) 


1.4343(20) 


3.0029(80) 


2.9447(26) 


0.90 


-1.5992(15) 


1.1742(22) 


2.820(13) 


2.7790(29) 


0.95 


-1.3887(13) 


0.8223(29) 


2.480(29) 


2.5639(36) 



TABLE XVII: Infrared-finite part of Qj for the unimproved 
action with kinetic term Hq only, as described in Ap- 
pendix IE II The gluon action is the Wilson action with 
A'^ = 10~^ and we use m — 2.8, n — 2. mean-field correc- 
tions are not included, the errors shown are statistical due to 
the VEGAS integration. 



FIG. 14: Heavy quark renormalization parameters for the 
simple, unimproved action with kinetic term Ho only. The 
gluon action is the Wilson action with A'^ — 10~® and we use 
m = 2.8, n = 2. All mean-field corrections are included and 
the results are infrared-finite. Note 6Zp is small due to small 
violation of reparametrization invariance. Also note SZm is 
large for this unimproved action. 



2. More improved case 



cellation described as in Appendix [Dl are (n = 2) 



(tadpole) 



(tadpole) 



.r(i-^ 



8 7TO 



3 - 2t>2 + ^4 



2 

^^(tadpolo) _ j=2(tadpolc) _ ^(2) ^ 2 I " 



16 

-9 + 3v^ 3 - 2v^ + 



47TO 



167774"^ 



(E2) 



The renormalization parameters of the heavy quark ac- 
tion (including mean-field corrections) are plotted in 

Fig. [131 For the one-loop coefficient of uq we use Mq = 
0.9735 [5l| and, as for the full action, we use cancellation 
of U^Uji^ described in the main text. 

We have also computed the renormalization parame- 
ters for the action discussed in [l9, [33| , 



iJo = -IV ■ A 
5H = 0, 



± A(2) - (-^ . A^)^ 
27TO 



where a different (less local) discretization of the operator 
{v ■ D)'^ is used. We used exactly the same simulation 
parameters as given there, m = 2.0, n = 2 and found 
agreement with their results for within statistical er- 
rors. 



We now consider a more improved action, including the 
spin-dependent cr ■ B'-term and the spatial and temporal 
lattice spacing improvement: 



SH = 



-iv ■ A 



± 



27TO 



(2) 



where SHcovi- is the same as in (1551) . We use the 
Symanzik-improved gluon action so that the Landau 
gauge mean link one-loop coefficient is Uq — 0.750 56J. 

Again, we considered m = 2.8 and n — 2, and as for the 
full 0{l/m'^,vf^i) action discussed in the main text the 
gluon mass was taken to be A^ = 10^®. The results for 
the rij obtained from the VEGAS integration are given in 
Table lXVlIII and we show the renormalization parameters 
(including mean-field corrections) as a function of the 
frame velocity in Fig. 1151 



APPENDIX F: POLES OF THE IMPROVED 
GLUON PROPAGATOR 

As the heavy quark action contains only first order 
time derivatives finding the poles in the propagator is 
trivial. It is also straightforward to find the poles of the 
simple, unimproved Wilson gluon propagator. However 
this is not the case for the Symanzik-improved gluon ac- 
tion. In this section we analyze the position of poles 
in the Symanzik-improved gluon propagator described in 
Ref. [ill. 
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2.1985(26) 
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-1.5932(14) 


1.3716(15) 


2.0890(29) 


1.9160(18) 
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-1.4613(13) 


1.2555(15) 


1.9695(35) 


1.7893(18) 
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-1.3094(12) 


1.1314(15) 


1.8435(46) 


1.6597(18) 
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FIG. 16: Poles of the Symanzik-improved propagator in the 
complex z plane. 



TABLE XVIII: Infrared-finite part of for the C»(l/m) 
action with chromomagnetic term, as described in Ap- 
pendix |E|2] The gluon action is the Symanzik-improved ac- 
tion with — 10~^ and we use m = 2.8, n = 2. Mean-field 
corrections are not included, the errors shown are statistical 
from the VEGAS integration. 
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frame velocity v 

FIG. 15: Heavy quark renormalization parameters for the 
0{l/m) action with chromomagnetic interaction term SH. 
The gluon action is the Symanzik-improved action with = 
10~® and we use m = 2.8, n = 2. All mean-field corrections 
are included. Note improvement drastically reduces 5Zm- 

We restrict our discussion to Feynman gauge where the 
gluon two-point function is given by 




+ (1 Qfiu)k^k^ 



with = 1 + ^{kl + kl) and fc^ = 2 sin(A:^/2). 

To find the poles of the propagator, first we compute 
the determinant of this matrix which is a polynomial in 



and Lo = kg. For a given three-momentum kj g [— tt, tt] 
the zeros of this expression in the z = e**^" plane can 
be obtained by solving det M (uj) = and then using 
uj = 2 — z — 1/z. It turns out that the determinant can 

be factored as det M(w) (w + fc + A^) det M{cj), with 

k ~ 4^^^j^ sin^(fcj/2), so that one solution coincides 
with the root of the naive propagator. Numerically, for 
small a^X^ also one of the solutions of detM(a;) = is 
very close to the naive solution. Note that the solutions 
come in pairs, (2;+,z_) with z^Z- = 1, so one of them 
lies inside the unit circle and the other outside. 

For a given spatial momentum there are 14 solutions. 
In Fig. [in] these are plotted in the complex z-planc for 
1000 randomly chosen kj. For the gluon mass a value of 
= 10^^ was chosen. 

To compare the poles in the improved propagator to 
the naive poles, their absolute value is computed and it 
is compared to that of the naive poles given by 

^£„aive) ^ l(^2 + fc' + A2 (F2) 

±V^(fc' + A2)(fe' + A2+4)) . 

In Fig. [T7] these absolute values are plotted for the same 
random three momenta. As can be seen from this plot the 
absolute value of an improved pole is either larger than 
^(naivc) smaller than 2:^"^*'™' but it never lies between 
these values. 

We performed a similar analysis for the propagator in 
Coulomb gauge and find that also in this case the poles 
of the Symanzik-improved propagator always lie outside 
the band defined by ^f"^'™^ < |z| < z^'^'™^ 

Hence, it is legitimate to use the position of the naive 
poles when deforming the integration contour in the de- 
termination of the heavy quark renormalization parame- 
ters. 
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FIG. 17: Absolute value of poles in the naive and 
Symanzik- improved gluon propagator as a function of — 
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